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ON 

VORTEX  SIMULATION  OF  TURBULENT  COMBUSTION 

(AFOSR  Grant  No.  89-0491) 

Principal  Investigator  Ahmed  F.  Ghoniem 

Department  of  Mechanical  Engineering 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  2139 


SUMMARY 

During  the  course  of  this  woik,  we  focused  on  the  extension  of  the  vortex  method  and  the 
transport  element  method  to  three-dimensional  flows  and  to  reacting  flow  with  Enite  heat  release, 
density  variation  and  volumetric  expansion.  Problems  explored  in  detail  include  the  formation  of 
streamwise  vorticity  in  uniform-density  and  variable  density  flow,  the  effect  of  this  vorticity  on 
the  rate  of  product  formation  and  the  effect  of  volumetric  expansion  due  to  heat  release  on  shear 
layer  growth  rate. 
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June  1990. 
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Progress  in  Astronautics  and  Aeronautics,  131,  ed.  by  Kulh,  Leyer,  Borisov  and 
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Entrainment  in  a  Shear  Layer,”  J.  Comput  Phys.,  97,  pp.  172-223,  Nov.  1991. 

6.  Ghoiuem,  A.F.  and  Heidarinejad,  G.,  “Effect  of  Two-dimensional  Shear  Layer  Dynamics 
on  Mixing  and  Combustion  at  Low  Heat  Release,”  Combust.  Sci.  and  Tech.,  Vol.  72,  pp. 
79-99,  March  1990. 

7.  Ghoiuem,  A.F.  and  Heidarinejad,  G.,  “Effect  of  Damkohler  Number  on  the  Reaction 
Zone  in  a  Reacting  Shear  Layer,”  Combust.  Flame,  83,  pp.  1-17, 1991. 

8.  Krishnan,  A.  and  Ghoniem,  A.F.,  “Simulation  of  the  Roll-up  and  Mixing  in  Rayleigh- 
Taylor  Flow  using  the  Vortex^’ransport  Element  Method,”  J.  Comput.  Phys.,  5^9, 1992, 
pp.  1-27. 

9.  Ghoiuem,  A.F.,  “Vortex  Simulation  of  Reacting  Shear  Flow,”  Chapter  10  in  Numerical 
Approaches  to  Combustion  Modelling,  Ed  by  E.  Oran  and  J.  Boris,  Progress  in 
Aeronautics  and  Astronautics,  Vol.  135, 1991,  pp.  305.348. 

10.  Ghoniem,  A.F.,  “Numerical  Simulation  of  Turbulent  Combustion  using  Vortex 
Methods,”  Fluid  Mechanical  Aspects  of  Combustion,  ed  by  M.  Onofri  and  A.  Resei, 
Pitman,  1991,  pp.  305-326. 

1 1.  Knio,  O.M.  and  Ghoniem,  A.F. ,  “Vortex  Simulation  of  a  Three-dimensional  Reacting 
Shear  Layer  with  Infinite-Rate  Kinetics,:  AIAA  Journal,  30,  January  1992. 

11.  Ghoniem,  A.F.,  Knio,  O.M.  and  Krishnan,  A.,  “Lagrangian  Simulation  of  the  Initial 
Stages  in  a  Reacting  Jet,:  The  23rd  Symposium  (International)  on  Combustion,  The 
Combustion  Institute,  Pittsburgh,  PA,  pp.  699-705, 1990. 

13.  Ghoniem,  A.F.,  Knio,  O.M.  and  Soteriou,  M.,  “Manipulation  of  the  Growth  Rate  of  a 
Variable-Density,  Spatially  Developing  Mixing  Layer  via  External  Modulation,” 
presented  at  the  29th  AIAA  Aerospace  Sciences  Meeting,  Reno,  NV,  January  1991, 

AIAA-91-0081. 


4 


I 


14.  Soteriou,  M.C.,  Knio,  O.M.,  Ghonietn,  A.F.  and  Cetegen,  B.M.,  “Simulation  of  Flow- 
Combustion  Interactions  in  a  Spatially  Developing  Mixing  Layer,”  presented  at  30th 
Aerospace  Sciences  Meeting  and  ExWbit,  January  6-9, 1992,  Reno,  Nevada,  AIAA-92- 
0081. 

15.  Ghoniem,  A.F.  and  Knio,  O.M.,  “The  Development  and  Application  of  the  Transport 
Element  Method  to  3D  Reacting  Shear  Layers,”  in  Vortex  Dynamics  and  Vortex 
Methods,  ed.  by  C.R.  Anderson  and  C.  Greengard.  Lectures  in  Applied  Mathematics, 

Vol.  28,  pp.  165-218, 1991. 

16.  Ghoniem,  A.F.,  Soteriou,  M.C.  and  Knio,  O.M.,  “Effect  of  Steady  and  Periodic  Strain  on 
Unsteady  Flamelet  Combustion,:  Proceedings  of  the  24th  Symposium  (International)  on 
Combustion,  July  5-10, 1992,  Sydney,  Ausmalia,  the  Combustion  Institute,  Pittsburgh, 
PA,  in  press. 

17.  Knio,  O.M.  and  Ghoniem,  “The  Three-din«nsional  Structure  of  Periodic  Vorticity  Layers 
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series  (by  invitation  only). 
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3.  1989,  National  Institute  of  Standards  and  Technology. 

4.  1989,  Ford  Motor  Company  Meeting  on  Flow  Visualization  (invited  speaker). 

5.  1990,  University  of  California,  Berkeley,  CA  (seminar). 

6.  1990,  California  Institute  of  Technology,  (seminar). 

7.  1990,  University  of  Southern  California,  (seminar). 

8.  1990,  NASA  Ames  Research  Center,  (seminar). 

9.  1990,  Sandia  National  Laboratories,  (seminar). 

10.  1990,  Lawrence  Livermore  National  Laboratories,  (seminar). 

1 1 .  1990,  Five-lecture  condensed  course,  Nagoya  University,  Japan 

12.  1990,  Institute  of  Computational  Fluid  Dynamics,  Tokyo,  Japan,  one-day  professional 
course. 

13.  1990,  Honda  Motor  Co.  Wako,  Japan,  (lecture). 

14.  1990,  Toyota  Motor  Co.,  Nagoya,  Japan,  (lecture). 

15.  1990,  Nissan  Motor  Co.,  Yokosaka,  Japan,  Qecture). 

16.  1990,  Mitsubishi  Heavy  Industries,  Nagoya,  Japan.(talk  and  discussion). 

17.  1990,  Mazda  Motor  Co.,  Yokohama,  Japan,  (talk  and  discussion). 

18.  1990,  American  Math  Society  Meeting  on  Vortex  Dynamics  (by  invitation  only),  Seattle. 

19.  1990,  SIAM  National  Meeting  (invited  talk),  Chicago. 

20.  1990,  USSR  academy  of  Sciences,  Novosibirsk,  USSR. 

21.  1990,  USSR  academy  of  Sciences,  Moscow,  USSR. 

22.  1991,  Gas  Research  Institute,  Chicago,  IL. 

23.  1991,  National  Center  for  Supercomputer  Applications,  Oiampaign,  IL. 

24.  1991 ,  ASME  Local  Chapter,  Boston,  MA. 
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25.  1991,  Northeastern  University,  Boston,  MA. 

26.  1991,  Ch.  Michelson  Institute,  Bergen,  Norway. 

27.  1992,  Egyptian  Environmental  Affairs,  (by  invitation)  Cairo,  Egypt. 

28.  1992,  Carnegie  Mellon  University  (seminar). 

29.  1992,  University  of  Massachusetts,  Amherst  (seminar). 

30.  1992,  Institute  of  Advanced  Studies  Workshop  on  Statistical  Physics  and  Fluid 
Mechanics,  Princeton,  NJ  (invited  lecture). 

31.  1992,  University  of  Tennessee  Space  Institute  Summer  Course  on  High  Angle  of  Attack 
and  Unsteady  Flow  Phenomena  (2  lectures). 
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TECHNICAL  REPORT 


The  transport  element  method  has  been  extended  to  three  dimensions  to  study  the 
evolution  of  scalar  fields  in  a  flow  with  high  vorticity  concentration.  The  numerical  scheme  is 
based  on  following  Lagrangian  computational  elements  employed  in  the  transport  of  vorticity 
and  local  scalar  gradients,  and  is  a  direct  generalization  of  the  three-dimensional  vortex  element 
method.  The  numerical  algorithms  required  to  implement  this  scheme  have  been  developed. 
Problems  associated  with  severe  distortion  of  the  flow  map  due  to  the  growth  of  perturbations 
were  shown  to  cause  difficulties  including  loss  of  numerical  accuracy  and  resolution.  Means  to 
overcome  these  problems  were  developed  and  were  shown  to  yield  accurate  solutions.  Two  grid- 
based  and  two  grid-free  methods  for  the  computation  of  vorticity  stretching  were  implemented, 
and  the  accuracy  of  the  methods  was  investigated  in  light  of  numerical  results.  Solutions  reveal 
the  need  for  a  careful  treatment  of  the  discrete  form  of  the  vorticity  transport  equation.  The 
methods  were  applied  to  study  the  evolution  of  an  initially  two-dimensional  shear  layer, 
perturbed  in  the  streamwise  and  spanwise  direction,  with  attention  focused  on  the  role  of  the 
spanwise  instability  and  how  it  enhances  the  rate  of  scalar  entrainment  into  the  large-scale 
structures  which  form  as  the  streamwise  instability  develops.  Two  mechanisms,  associated  with 
the  onset  of  three-dimensional  instability,  are  found  to  be  responsible  for  this  enhancement: 
vorticity  intensification  within  the  large  eddy  core  due  to  spanwise  stretching  which  delays  its 
collapse;  and  generation  of  transverse  entrainment  currents  towards  the  eddy  core  due  to  the 
formation  of  streamwise  vortex  structures  within  the  core  and  along  the  braids  between 
neighboring  cores.  Preferential  entrainment  is  detected  along  the  spanwise  direction  due  to  the 
streamwise  vorticity.  Details  of  this  work  are  presented  in  the  paper  on,  “Three-dimensional 
Vortex  Simulation  of  Rollup  and  Entrainment  in  a  Shear  Layer,”  by  Omar  M.  Knio  and  Ahmed 
F.  Ghoniem.(i) 

Numerical  simulations  of  a  three-dimensional  temporally-growing  shear  layer  were  then 
obtained  at  high  Reynolds  number  for  a  variable  density  flow  using  the  transport  element 
method.  Attention  was  focused  on  the  effect  of  initial  vorticity  and  density  distributions  on  the 
interaction  between  the  instability  modes  which  lead  to  the  generation  and  intensification  of 
streamwise  vorticity.  Results  showed  that  the  three-dimensional  instabilities  evolve  following 
the  formation  of  concentrated  spanwise  vorticity  cores.  The  deformation  of  each  core  along  its 
span  resembles  the  amplification  of  the  translative  instability.  The  generation  of  vortex  rods 
which  wrap  around  individual  cores  while  stretching  between  neighboring  cores,  suggested  a 
mode  similar  to  the  Corcos  instability.  The  instability  modes  leading  to  the  formation  of  both 


(1)  Journal  of  Computational  Physics,  Vol  97,  No.  1,  November  1991. 
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structures,  energized  by  the  extensional  strain  generated  by  the  cores,  grow  simultaneously.  A 
similar  series  of  events  occurred  in  variable-density  shear  layers  and  in  shear  layers  which  start 
with  an  asymmetric  vorticity  distribution.  Baroclinic  vorticity  generation  in  the  variable-density 
layer  lead  to  the  formation  of  asymmetric  cores  whose  volumetric  composition  is  biased  towards 
the  lighter  fluid.  The  structures  are  propelled,  by  their  asymmetric  vorticity  distribution,  in  the 
direction  of  the  heavier  stream  while  their  eccentric  spinning  forces  an  uneven  stretching  of  the 
vortex  rods.  The  origin  of  the  asymmetry  was  established  by  comparing  the  results  of  variable 
density  shear  layers  with  the  results  of  a  shear  layer  with  an  initially  asymmetric  vorticity 
distribution  in  a  uniform-density  flow.  The  strong  late-stage  asymmetry  exhibited  by  the  former 
is  not  observed  in  the  latter.  Thus,  baroclinic  vorticity  generation  is  responsible  for  the  observed 
asymmetry.  We  also  find  that  initially  asymmetric  vorticity  distribution  does  not,  as  suggested 
before,  lead  to  asymmetric  spacing  between  the  streamwise  rods.  It  is  concluded  that  the 
experimentally  observed  asymmetric  spacing  must  arise  after  pairing.  Details  of  this  study  are 
presented  in  this  report  as  Appendix  I.  (2) 

The  three-dimensional  transport  element  method  was  then  extended  to  solve  the 
conservation  equations  for  reacting  flow.  The  numerical  scheme  maintained  its  adaptive, 
Lagrangian  nature  in  which  computational  effort  is  concentrated  in  zones  of  finite  vorticity  and 
chemical  reaction.  The  method  utilized  an  accurate  discretization  of  vorticity  and  species 
concentrations  among  a  number  of  spherically-symmetric  reacting  transport  element  and  the 
advection  of  these  elements  along  particle  trajectories.  We  used  the  low  Mach  number 
approximation  of  combustion  in  open  domains  and  restricted  our  attention  to  the  case  of 
diffusion  flames  with  no  heat  release.  A  single-step,  second-order,  infinite-rate  kinetics  chemical 
reaction  model  was  employed.  The  scheme  was  applied  to  study  the  effect  of  flow-induced 
instabilities  on  the  reaction  zone  and  product  distribution  in  a  temporal  shear  layer.  Results  were 
obtained  in  the  high  Peclet  number  regime  for  a  wide  lange  of  Damkohler  numbers.  Changes  in 
the  reaction  field  were  related  to  either  the  entrainment  or  the  strain  field  associated  with  the 
saturation  of  the  instabilities.  With  increasing  Damkohler  number,  the  structure  of  the  reaction 
region  changed  from  a  distributed  zone  embedded  within  spanwise  and  streamwise  vortices  to  a 
thin  sheet  surrounding  their  cores.  However,  the  product  concentration  exhibited  strong 
similarity  to  the  vorticity,  falling  rapidly  in  regions  where  the  vorticitv  magnitude  was  small. 
Variation  of  the  Peclet  number  yields  minor  changes  in  the  product  u. aibution  and  in  the 
reaction  zone  structure  but  strongly  affects  product  formation  rates.  Details  of  this  study  are 
printed  in  this  report  as  Appendix  II. 


The  paper  has  been  accepted  for  publication  in  the  Journal  of  IHuid  Mechanics. 
(^1  The  paper  has  appeared  in  the  AIAA  Journal,  Vol  30,  pp.  105- 1 16,  January  1992. 
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Finally,  the  transport  element  method  has  been  extended  to  simulate  a  spatially 
developing,  reacting  shear  layer  with  unmixed  reactants  and  finite  heat  release.  In  the  case  of  a 
reacting  flow,  vorticity  changes  due  to  density  variation,  scalar  gradients  change  due  to  the 
chemical  reaction,  and  volumetric  expansion  adds  an  expansion  field.  Solutions  were  obtained 
for  a  forced  shear  layer  at  different  Damkohler  numbers  and  enthalpy  of  reaction  to  study  the 
effect  of  combustion  heat  release  rate  on  the  development  of  the  large  scale  structures.  Forcing 
was  used  to  ensure  roll-up  within  the  computational  domain.  We  showed  that  heat  release 
enlarges  the  size  of  the  fundamental  eddies,  stretching  their  streamwise  dimension  and  slightly 
reducing  their  cross  stream  dimension,  while  their  overall  size  remains  almost  the  same.  Along 
with  forcing  at  the  fundamental  and  subharmonic  frequencies,  heat  release  increases/decreases 
the  size  of  the  in-phase/out-of-phase  eddies.  The  non-uniform  acceleration  of  the  eddies  in  the 
streamwise  direction  causes  their  relative  locations  to  deviate  from  that  of  a  uniform-density 
layer  and  thus  modifies  the  pairing  process  into  a  tearing/gulping  process.  Results  also  show  that 
the  enthalpy  of  reaction  is  more  important  than  the  reaction  frequency  factor  in  affecting  the  flow 
dynamics.  For  the  same  parameters,  the  variable-density  layer  grows  slower  than  its  uniform- 
density  counterpart.  Details  of  this  work  are  presented  in  Appendix  in  in  this  report. 


Presented  at  the  30th  Aerospaace  Sciences  Meeting  and  Exhibit,  Jan  6-9, 1992,  Reno,  NV,  AIAA-92-0081. 
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APPENDIX  I 


THE  THREE-DIMENSIONAL  STRUCTURE  OF  PERIODIC  VORTICITY  LAYERS 
UNDER  NON-SYMMETRIC  CONDITIONS 


Omar  M.  Knio^  and  Ahmed  F.  Ghoniem 
Dq)aitment  of  Mechanical  Engineering 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 


ABSTRACT 

Numerical  simulations  of  a  three-dimensional  ten^orally-growing  shear  layer  are  obtained 
at  high  Reynolds  number  and  zero  Fronde  nunriier  uring  a  vortex  schone  modifi^  for  a  variable 
density  flow.  Attention  is  focused  on  the  effect  of  initial  vorticity  and  density  distributions  on  the 
interaction  betweoi  instalnlity  modes  which  lead  to  die  generation  and  intensificati<m  of  streamwise 
vorticity.  Results  show  that  the  three-dimoisional  instabilities  evolve  following  the  formation  of 
concentrated  spanwise  vorticity  cores.  The  defamation  of  each  core  along  its  span  resembles  the 
amplification  of  the  translative  instability.  The  generation  of  vortex  rods,  wMch  wrap  around 
individual  cores  while  stretching  between  neighbe^g  cores,  suggest  a  mode  similar  to  Ae  Coicos 
instability.  The  instability  modes  leading  to  die  foimadon  of  both  structures,  ena*gized  by  the 
extension  strain  generated  by  the  cores,  grow  rimultaneously.  A  similar  series  of  events  occurs 
in  variable-density  shear  layers  and  in  shear  layers  which  start  with  an  asymmetric  vorticity 
distribution.  Baroclinic  vorticity  generatiem  in  die  variable-density  layer  leads  to  die  formation  of 
asymmetric  cores  whose  volumetric  contrition  is  Inased  towards  the  lighter  fluid.  The  structures 
are  propelled,  by  dieir  asymmetric  vorticity  distribotion,  in  the  direction  of  the  heavier  stream  while 
their  eccentric  spinning  forces  an  uneven  stretching  of  the  vortex  rods.  The  origin  of  the 
asymmetry  is  establish^  by  conqiaring  these  with  the  results  of  a  shear  layer  with  an  initially 
asymmetric  vorticity  distribution  in  a  uniform-density  flow.  The  strong  late-stage  asymmet^ 
exhibited  by  the  former  is  not  observed  in  the  latter.  Thus,  baroclinic  vorticity  generation  is 
responsible  for  the  observed  asymmetry.  We  also  find  that  initially  asymmetric  vorticity 
distribution  does  not,  as  suggested  before,  lead  to  asymmetric  spacing  between  the  streamwise 
ro^.  It  is  concluded  that  the  experimmtadly  observed  asymmetric  spacing  must  arise  after  pairing. 
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1.  INTRODUCTION 


The  fomiation  of  large  vortical  structures  has  long  been  observed  in  free  shear  layers  at 
high  Reynolds  numbers  (Crow  &  Champagne,  1970;  Brown  &  Roshko,  1974).  Analysis  of 
experimental  results  shows  that  the  evolution  of  these  structures  and  their  mutual  interactions, 
governed  essentially  by  inviscid  flow  dynamics,  play  an  important  role  in  die  growth  of  the  layer, 
the  distribution  of  turbulent  statistics,  scalar  transport  and  mixing  (Winant  &  Browand  (1974);  Ho 
&  Huerre,  1984;  Dimotakis,  1989;  Hussain,  1986).  Mthout  external  forcing,  the  early  stages  of 
the  layer  ate  dominated  by  2D  modem,  where  spanwise  vortex  structures  are  formed,  followed  by  a 
transidon  to  3D  modon.  Within  and  after  the  transidon  region,  the  spanwise  vortices  are 
deformed,  and  secondary  streamwise  vortices  are  generated  (Konrad,  1976;  Breidenthal,  1981; 
Jimenez,  1983;  Bernal,  1981).  Experimental  investigations  yielded  striking  visualizations  of  the 
streamwise  structures,  and  how  they  modify  scalar  transport  (Jimenez,  Cogollos  &  Bernal,  1985; 
Bernal  &  Roshko,  1986).  In  recent  work,  Lasheras,  Cho  &  Maxworthy  (1986)  and  Lasheras  & 
Choi  (1988)  exaiitined  the  possibility  of  manipulating  the  location  of  the  "transition"  re^von  by  3D 
forcing  (see  also  Breidenthal,  1980),  thereby  eti^hasizing  die  practical  aspect  of  such  studies. 

The  3D  response  of  votticity  layers  is  con^lex  and  determining  the  origin  and  shape  of  die 
secondary  structures  poses  considerable  difficult^  to  analytical  studies.  So  far,  a  limited  number 
of  dieotetical  investigations  have  dealt  with  this  problon.  Pkntdiumbert  and  Widnall  (1982)  used 
a  periodic  airay  of  Stuart  vortices  to  represrat  the  ^anwise  eddies  formed  by  the  roll>up  of  the 
Kelvin-Helmholtz  instability.  The  linear  stability  analysis  of  this  configuration  revealed  the 
presence  of  a  "translative"  instability,  which  was  thm  proposed  as  a  possible  mechanism  leading  to 
the  formation  of  the  observed  secondary  vortices  in  shear  layers.  Anodier  mechanism,  discovoed 
by  Coicos  &  Un  (1984)  in  dieir  study  of  the  staUlity  of  the  layer  linearizing  3D  perturbations 
around  die  evolving  2D  flow  (Lin  &  Coicos,  1984),  was  also  suggested.  They  showed  that,  for 
sufficiendy  low  difiiision  (Neu,  1984),  the  strained  streamwise  vorticity  is  unstable,  and  the 
instability  causes  the  redistribution  of  the  latto’  into  round,  concentrated  vortex  rods.  These 
streamwise  rods  lead  to  the  generation  of  "mushroom"  structures  similar  to  those  experimentaUy 
observed  (Bernal  &  Roshko,  1986;  Lasheras,  Cho  &  Maxwortfay,  1986;  Lasheras  &  Choi,  1988). 
Manifestation  of  bodi  instabiliQr  mechanisms  has  been  npoittd  in  experimental  studies,  however, 
the  broad-band  nature  of  die  3D  modes  predicted  in  both  theories  complicates  the  task  of  verifying 
thdr  validity  or  determining  flow  conditions  under  which  one  mode  dominaies  die  odier. 

The  3D  motion  of  shear  layers  has  also  been  the  subject  of  numerical  investigations. 
Ashurst  &  Meiburg  (1988)  used  a  vortex  filament  scheme  to  conqiute  the  development  of  a 
temporal  shear  layer  at  high  Reynolds  number.  The  shear  layer  was  modeled  by  a  single 
desingularized  vortex  sheet  <x  two  vortex  sheets  of  opposite  sign.  Results  oi  both  models  showed 
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evidence  of  both  the  translative  and  the  Corcos  instability  modes.  When  two  sheets  of  opposite 
signs  were  used,  an  asymmetric  distribution  of  streamwise  vortices  (Lasheras  &  Choi,  1988),  was 
apparently  reached  through  a  non-linear  interaction  between  two  counter-rotating  streamwise 
vortex  rods,  each  originating  in  a  distinct  vorticity  layer.  While  a  similar  asymmetric  distribution 
of  vortices  was  observed  experimentally  in  shear  layers,  the  initial  vorticity  profile  used  in  this 
simulation  was  more  representative  of  that  of  a  wake.  The  difference  between  the  stability  and 
long-time  behavior  of  wakes  and  shear  layers  could  be  seen  by  cong>aring  the  results  of  Ashurst  & 
Meiburg  and  of  Grinstein  et  al.  (1989),  who  used  a  finite-difference  scheme  to  simulate  the 
evolution  of  a  spatially-developing  shear  layer  at  moderate  Reynolds  number.  Results  of  the  latter 
indicate  that  the  asymmetric  distribution  of  streamwise  vortices  develops  as  a  result  of  the  merging 
of  pairs  of  streamwise  rods  of  the  same  sign  of  circulation,  and  that  this  interaction  only  occurs 
after  pairing  between  neighboring  spanwise  eddies. 

The  "sequential"  nature  of  the  growth  of  several  forms  of  2D  and  3D  instabilities  was 
investigated  in  the  spectral  calculations  of  Metcalfe  et  al.  (1987).  They  considered  temporal 
vorticity  layers  at  low  Reynolds  number,  and  performed  a  detailed  study  of  the  energy  content  of 
the  distinct  modes.  Their  results  show  that,  for  small  perturbations,  2D  Kelvin-Helmholtz  waves 
grow  first  During  their  growth,  3D  activity  is  suppressed  and  the  layer  rruuntains  a  2D  charactCT. 
Following  the  non-linear  growth  of  the  Kelvin-Helmholtz  waves  and  the  formation  of  ("primary") 
spanwise  eddies,  3D  perturbations  are  amplified.  Saturation  of  tire  Kelvin-Helmholtz  and  3D 
instabilities  is  reached  soon  after  the  flattening  of  the  (primary)  spanwise  cores  and  the 
"maturation"  of  3D  modes.  They  also  show  that  pairing  of  spanwise  cores,  during  which  3D 
activity  is  suppressed,  is  necessary  for  further  growtii  of  the  layer.  Following  pairing,  the  growth 
of  3D  modes  is  resumed.  An:q>lification  of  3D  disturbaiKes  is  thus  restricted  to  windows 
separating  periods  of  otherwise  2D  growth. 

The  dependence  of  tire  response  of  the  shear  layer  on  initial  conditions  and  forcing  levels 
was  investigated  by  Inoue  (1989),  who  performed  a  vortex  filament  simulation  of  a  spatially- 
developing  layo*  modeled  by  a  single  vortex  sheet  His  results  show  that  the  3D  transition  strongly 
depends  on  3D  forcing,  and  that  the  flowfield  terxls  towards  2D  behavior  once  this  forcing  is 
interrupted.  This  was  contradicted  by  the  numerical  experiments  of  Grinstein  et  al.  (1989)  which 
indicate  that  the  transition  to  3D  motion  persists  despite  the  absence  of  a  continuous  forcing 
function.  This  result  was  verified  in  the  simulations  of  Lowery  et  al.  (1987),  who  employed  a 
hybrid  finite  difference  -  q)ectral  method  to  track  the  evolution  of  a  passive  scalar  in  a  developing 
layer  at  low  Reynolds  number.  While  the  study  focused  on  asymmetric  entrainment  patterns  in  2D 
and  3D  (Dimotakis  1989, 1986),  it  also  emphasized  the  fact  that  distributions  of  ^anwise  and 
streamwise  vorticity  are  weakly  dependent  on  the  strength  and  sh^  of  the  forcing  function.  As 
indicated  below,  vortex  simulations,  similar  to  those  of  Inoue  (1989)  and  Ashurst  &  Meiburg 
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(1988),  which  do  not  allow  the  number  of  con^utational  elements  to  increase  as  the  material 
elements  are  stretched,  may  lose  accuracy  and  contaminate  the  computations  with  numerical 
di^sion  errors. 

Despite  these  efforts,  a  clear  and  unified  interpretation  of  the  role  of  the  various  modes  in 
the  formation  of  vortical  and  scalar  structures,  a  crucial  step  towards  a  bettor  understanding  of  3D 
transition,  has  not  been  reached.  Moreover,  sevoal  important  issues,  summarized  next,  cast  some 
doubt  on  die  conclusions  of  these  simulations.  Using  de^gulaiized  vortex  sheets  to  model  shear 
layers,  as  in  Ashurst  &  Meiburg  (1989),  may  lead  to  spurious  results  since,  as  shown  by  the 
detailed  numerical  study  of  Knio  &  Ghoniem  (1991),  the  properties  of  the  3D  modes  of  a  vortex 
structure  are  strongly  dependent  on  the  vorticity  distribution  within  the  cross  section  of  the 
structure.  The  use  of  vortex  filaments  to  resolve  vorticity  within  the  shear  layer  is  not 
recommended  since  the  cotiesponding  schemes,  although  accurate  for  short  times,  do  not  maintain 
their  accuracy  as  vortex  elements  tend  to  move  apart  due  to  stretching  in  die  direction  normal  to  that 
of  die  main  flow.  This  violates  an  in^rtant  accuracy  condition  in  vortex  methods,  namely  diat  the 
cores  of  nei^boring  vortex  elements  must,  at  all  times,  overiiq). 

Another  source  of  inaccuracy  in  the  simulations  presented  in  Ashurst  &  Meiburg  (1988)  is 
the  enforcement  of  the  periodic  boundary  conditions.  A  small  number  of  images  was  used  on 
either  side  of  the  computational  domain,  which  was  not  enough  to  cipture  the  correct  value  of  the 
free  stream  velocity.  On  the  other  hand,  the  conclusions  of  the  spectral  simulations  of  Metcalfe  et 
al  (1987),  particularly  those  concerning  the  onset  of  and  the  interaction  between  the  2D  and  3D 
inodes  described  above,  may  be  dqiendent  on  some  diffusive-convective  balance  achieved  only  at 
low  Reynolds  number.  This  balance,  which  is  a  function  of  the  Reynolds  numbo-,  can  lead  to  an 
early  saturation  of  the  instability.  We  believe  that  a  hi^  Reynolds  number  siimilation  is  necessary 
to  determine  whether  the  intoractitm  between  die  different  modes  of  the  3D  instabiliQ'  is,  as  widely 
suspected,  an  essentially  inviscid  iiKchanisia 

The  large  number  of  mechanisms  governing  the  evolution  of  shear  layers,  and  the 
complexity  of  die  resulting  vortical  and  scalar  stmOures  underscore  the  need  for  accurate  numerical 
methods  which  can  carefully  treat  the  the  vorticity  transport  equation.  Successful  implementation 
of  numerical  schemes  depends  on  the  inoper  account  of  the  vorticity  stretching  term,  and,  if 
present,  vorticity  source  terms.  Another  crucial  ingredient  lies  in  the  ability  of  the  method  to 
accommodate  the  large  strain  associated  widi  high  concentratitms  of  vorticity.  The  latter  have  been 
shown  to  result  in  the  deterioration  of  the  discretization  accuracy  in  both  Euloian  computations, 
through  die  creation  of  small  scale  structures  which  may  not  be  well  represented  on  a  grid  of  fixed 
mesh  size  and  the  accumulation  of  numerical  diffusion  which  may  dissqiato  these  structures  at  their 
early  stages,  and  Lagrangian  computations  where  zones  of  high  strain  may  be  depleted  of 
computational  elements.  Finally,  questions  regarding  the  dynamic  effect  of  density  gradients,  in 
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the  absence  of  gravity,  which  impact  mixing  of  gases  at  different  molecular  weight  and/or  density, 
have  not  been  tackled  before. 

In  this  work,  an  adaptive,  Lagrangian  numerical  scheme  is  used  in  the  simulation  of 
vorticity  layers.  Two  main  ingredients  are  incorporated  in  the  construction  of  the  scheme,  which 
was  analyzed  in  our  previous  effort  (Knio  &  Ghoniem,  1991).  The  first  relates  to  its  adaptive 
nature,  a  feature  which  avoids  the  loss  of  spatial  resolution  and  allows  the  accommodation  of  high 
strain  rates  by  increasing  the  number  of  computational  elements  as  the  flow  evolves.  The  second, 
a  properQ^  found  in  most  Lagrangian  methods,  consists  of  ensuring  a  minimal  effect  of  numerical 
diffusion  which  may  lead  to  excessive  smearing  of  the  vorticity.  We  also  apply  a  vortex  scheme 
extended  to  variable-density  flows  to  analyze  the  dynamic  effect  of  finite  density  gradients  on  the 
evolution  of  the  shear  layer. 

Cmnputed  results  are  used  to  accurately  portray  die  severe  deformation  of  die  flow  map  and 
the  evolution  of  the  flow  vorticity.  We  focus  on  the  relationship  between  the  deformation  of 
material  surfaces,  the  generation  and  intensification  of  vorticity  and  the  associated  scalar 
entrainment  patterns.  The  results  are  used  to  characterize  die  3D  instabilities  of  the  vorticity  layer. 
Instability  modes  leading  to  the  generation  of  stieamwise  vortex  rods  joining  neighboring  eddies 
are  identified  and  distinguished  from  those  affecting  die  v(»rtex  cores.  A  generalized  perspective  of 
the  latter  is  given;  the  similarity  between  vorticity  patterns  found  in  the  late  stages  of  spanwise 
vortex  cores  in  the  shear  layer  and  those  observed  in  the  develt^rnmat  of  vortex  rings  is  discussed. 
Finally,  we  investigate  some  of  the  mechanisms  leading  to  the  onset  of  asymmetry  in  an  otherwise 
symmetric  flow.  In  particular,  a  variable-density  layer  is  contrasted  with  a  uniform-density 
asymmetric  layer,  in  order  to  study  the  roles  of  density  variation  and  asymmetric  strain  field  on  the 
development  of  the  vorticity  field. 

The  numerical  scheme  used  in  die  computations  is  summarized  in  Section  2.  In  Section  3, 
we  review  the  evolution  of  a  uniform-density  symmetric  layer,  then  present  results  of  variable- 
densi^  and  asymmetric  vorticity  layer  coiiqiutations.  The  results  are  further  discussed  in  Section 
4;  concluding  mnarks  are  given  in  Section  5. 

2.  FORMUIJ^TION  AND  NUMERICAL  SCHEME 

2.1.  FORMULATION  AND  GOVERNING  EQUATIONS 

We  start  with  the  incompressible,  isentropk,  variable-density  form  of  the  governing 
equations  in  the  low  Mach  number  limit  (Rehm  &  Baum,  1978;  Majda  &  Sethian  1987,  Ghoniem 
&  Krishnan,  1989).  By  choosing  an  appropriate  combination  of  characteristic  length,  time  and 
velocity  scales  as  normalizing  parameters,  these  equations  are  written  as: 
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(1) 


The  vorticity  associated  with  a  material  particle  changes  due  to  tilting  and  stretching  under  the 
action  of  the  strain,  Vu,  and  due  to  the  particle  acceleration  in  a  non-uniform  density  field. 
Equation  (S)  is  preferred  over  Eq.  (4)  because  the  baioclimc  source  term  is  written  in  tenns  of  the 
Idnanadcs  of  the  flowfield  ratter  than  being  treated  as  a  dynamic  effect  associated  with  pressure 
forces.  Buoyancy  effects  have  been  neglected  in  the  vorticity  tranqxat  equation  since  we  intend  to 
focus  on  high-speed  flows  in  which  fluid  acceleration  is  much  larger  than  gravitational 
acceleration.  Under  these  conditions,  gravity  effects  are  negligible  whenever  the  Richardson 
number, 

po  (6) 
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written  in  terms  of  the  gravity  constant,  g,  and  characteristic  density,  po ,  length,  Ai,  density 
difference,  Ap,  and  velocity  difference,  U,  is  small  (Koop  &  Browand,  1979). 

The  presence  of  the  baroclinic  source  tom  requires  the  accurate  estimation  of  the  density 
gradient  To  this  end,  we  derive  a  transport  equation  for  the  scalar  gradient  (the  density  in  this 
case),  by  taking  the  gradient  of  Eq.  (1)  to  obtain: 

Dg  „ 

-S-st .  o-vu  -  ex  CD 

Dt  ^  ^  ^  (7) 

where  g  =  Vp.  Thus,  while  density  remains  constant  along  a  material  path,  its  gradient  is  affected 
by  the  local  strain  and  the  voiticity.  Working  widi  Eq.  (7)  instead  of  Eq.  (1)  is  of  similar  nature  as 
using  the  voiticity  transport  equation  in  place  of  the  momoitum  equation.  Both  substitutions  are 
motivated  by  the  observation  that,  in  most  high  Reynolds  number  flows,  the  supports  of  the 
voiticity  and  scalar  gradient  are  small  subsets  of  die  supports  of  the  primitive  variables.  Thus, 
computational  effort  is  concentrated  into  smalla  regions  of  the  domain  of  study,  and  the  numerical 
differentiation  of  the  density  field  is  avoided,  thereby  mudmizing  a  loss  of  resolution. 

2.2.  NUMERICAL  SCHEME 

The  transport  element  method  is  used  to  confute  the  evolution  of  die  shear  layer.  The 
numerical  scheme  solves  the  time-dependent,  inviscid,  incompressible,  vorticity,  scalar  and  scalar 
gradient  transport  equations  given  above.  Variants  of  the  numerical  scheme  which  accommodate 
gravity,  compressibility,  chemical  reactions  and  diffusion  have  been  extensively  used  in  2D 
(Ghoniem  &  Kiishnan,  1988),  and  in  a  limited  number  of  applications  in  three-dimensions  (Knio 
&  Ghoniem,  1992),  but  will  not  be  required  in  this  study. 

The  numerical  method  is  the  product  of  the  combination  of  a  series  of  refinements  of  3D 
vortex  methods  with  the  scalar  tranq;K)rt  techniques  developed  in  the  2D  tian^rt  element  method 
(Ghoniem  et  al.,  1988).  It  is  based  on  the  discretization  of  the  vortici^  and  scalar  gradient  fields 
into  a  finite  numbor  of  Lagrangian  elements,  called  transport  elements.  Transport  elements  carry 
discrete  scalar,  votticity  and  scalar  gradient  values,  and  are  distributed  along  elemoitary  rectangular 
areas  which  are  used  to  divide  entire  material  surfaces.  The  Lagrangian  mesh  defines  the  location 
of  the  elements,  iKdiile  vector  quantities  are  not  restricted  to  lie  within  dementary  rectangular  areas. 
Discrete  quantities  are  snnootiied  in  a  small  spherical  neighborhood  of  the  centO’  of  die  element  A 
third  order  Gaussian  core  function,  yr;  s  3l(4nS^)  exp(-r^l&)^  is  used  as  smoothing  function,  and 
its  standard  deviation  5  is  used  as  characteristic  core  radius  (Knio  &  Ghoniem,  1991).  This  yields 
a  continuous  veraon  of  die  voiticity  field  which  induces  a  desingularized  velocity  field,  expressed 
in  terms  of  die  Biot-Savart  law  (BatchelOT,  1967).  The  material  surfaces  are  tracked  by  moving  the 
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vertices  of  the  elements  with  the  local  velocity  vector,  using  a  second-order,  predictor-corrector 
tune  integration  scheme. 

The  construction  of  the  transport  element  method  closely  mimics  conventional  vortex 
element  schemes.  Numerical  analysis  of  these  schemes  reveals  that  the  smoothing  functions 
control  the  order  of  the  spatial  resolution  and  that  strong  oveiiiq>  among  die  cores  of  neighboring 
elements  is  required  to  guarantee  convergence  (Beale  &  Majda,  1982a,b;  Beale,  1986).  While  the 
third-order  Gaussian  functions  have  been  shown  to  yield  second-order  schemes  (Leonard,  1985; 
Beale  &  Majda,  1985),  it  is  the  overlap  condition  that  first  motivated  the  modification  of  the  3D 
VOTtex  method,  in  which  vortex  elements  are  distributed  along  vortex  tubes.  In  vortex  methods 
and  vortex  filament  methods,  vortex  elements  are  redistributed  along  vortex  tubes  whatever  the 
strain  causes  the  separation  distance  between  neighboring  elements  to  exceed  the  core  radius. 
However,  the  overlap  condition  cannot  be  enforced  if  the  elements  (or  filaments)  are  strained  in  a 
directicMi  normal  to  die  local  vorticity  vector.  As  a  result,  deterioration  in  the  ^atial  resolution  may 
occur  at  long  times  as  large  discretizatitm  enrors  pollute  the  solution.  In  the  present  computations, 
a  scheme  of  local  mesh  refinement  which  subdivides  computational  elonents  along  two  directions 
of  strain  is  employed.  This  scheme,  which  is  described  below,  has  been  shown  to  yield  significant 
improvement  in  the  accuracy  of  the  confutations  over  vortex  element  computations  (Knio  ft 
Ghoniem,  1991). 

Another  advantage  of  the  transport  element  method  lies  in  die  fact  diat  tracking  material 
surfaces  greatly  siiif  lifies  the  task  of  integrating  the  equation  of  motion  of  the  density  gradient 
This  is  achieved  associatir<  vidi  each  tranfort  elonent  an  elementary  surface  area  SA.f(t)  and  a 
unit  normal  to  the  surface  area  at  the  center  of  the  element  Hi(t)  =  dAi(t)llSAi(t)l.  In  the 
confutations,  this  is  done  by  requiring  diat  material  surfaces  effectively  constitute  iso-scalar 
surfaces,  and  adopting  linear  interpolation  functions  to  describe  the  slufe  of  the  material  surface 
within  each  tranfort  element  We  take  advantage  of  Idnematical  relationships  which  relate  the 
evolution  of  the  gradient  of  a  non-diffusive  scalar  in  an  incompressible  fluid  at  a  material  point 
X(t),  g(t,x),  to  that  of  an  elementary  surface  area  coitered  around  x*  SA(t,x),  initially  having  same 
sense  and  direction  as  g.  This  relationship,  which  constitutes  the  analogue  of  die  Helmholtz 
vOTticity  theorem,  may  be  expressed  as:  g{t,x)  -  aiSA(t,xh  ^iriiere  a  s  lg(0,x)l/lSA(0,x)l‘  In  view 
of  the  preceding,  we  avoid  integrating  Eq.  (7)  simply  by  following  die  evolution  die  elementary 
surface  areas. 

The  description  of  the  method  is  cotif  leted  by  specifying  a  tBchnk}oe  for  updating  the 
discrete  values  of  vorticity  associated  widi  the  tranftnt  elements.  Two  techniques  are  used  in  the 
confutations.  In  die  absence  of  density  variation,  direct  integration  of  the  vorticity  transport 
equation  is  avoided.  In  this  case,  substantial  computatitmal  savings  are  achieved  1^  taking 
advantage  of  the  Helmholtz  and  Kelvin  theorems  since  vorticity  lines,  identified  by  their 
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circulation,  evolve  as  material  lines.  In  this  variant  of  the  scheme,  vorticity  changes  according  to 
the  stretching  and  tilting  of  elementary  material  segmoits  lying  in  its  direction,  while  the  circulation 
associated  with  a  transport  element  remains  constant  This  procedure  avoids  the  evaluation  of  the 
velocity  gradient,  and  retains  the  conservation  of  the  volume  of  vorticity  property  possessed  by 
vortex  filament  methods  (Greengard,  1986). 

When  dealing  witii  a  variable  density  flows,  we  can  no  longer  iq>ply  the  Helmholtz  vorticity 
theorem.  Therefore,  we  are  not  able  to  avoid  the  integration  of  Eq.  (5),  and,  in  doing  so,  the 
evaluation  of  the  velocity  gradient  The  procedure  suggested  here  is  to  q>lit  that  task  in  two 
firactional  steps,  by  first  numerically  integrating: 


fia=o).Vu 

Dt 

and,  in  a  second  step,  integrating: 

vDm 

Dt  p  Dt 

This  procedure  is  similar  to  the  "viscous  splitting"  of  the  viscous  vorticity  transport  equation 
(Qiorin,  1973),  and  thus  may  be  termed  "baroclinic  splitting"  of  the  equation  of  motion.  In  the 
numerical  integration  of  Eq.  (8),  Vu  is  found  by  analytically  differentiating  the  desingularized 
Biot-Savart  law,  while  the  same  predictor-corrector  onployed  to  advance  the  computational  mesh 
is  used  to  perform  die  time  integratitm.  Equation  (9)  is  integrated  in  a  single  step,  estimating  the 
baroclinic  torque  from  the  knowledge  of  (i)  p  which  is  constant  for  each  element,  (ii)  Vp  which  is 
computed  according  to  the  defoimation  of  the  elements,  and  (iii)  Du/Dt  which  is  approximated  by  a 
first-order,  backward  finite  diffoence  in  time,  DulDt  s(u(t)’U(t-^))ldt. 

In  the  non-linear  evolution  of  the  flowfield,  a  strong  and  rapid  deformation  of  the 
Lagrangian  computatimial  mesh  is  experienced.  This  deformation  causes  the  depletion  of 
conq>utational  elements  in  some  regions  of  the  domain  where  the  separation  distance  between 
neighboring  elements  becomes  excessively  large.  In  this  work,  we  employ  a  local  mesh 
refinement  scheme  which  splits  a  transport  element  into  two  whenevCT  the  average  value  of 
opposing  sides  of  the  rectangles  exceeds  the  value  of  the  core  radius.  The  scheme  has  been 
described  in  detail  in  (Knio  &,  Ghoniem,  1991).  It  essentially  ensures  tiut  core  overlap  annong 
iwi^boring  elements  is  satisfied  and  amounts  to  redistributing  the  vorticity  and  scalar  fields  into  a 
larger  number  of  elements  due  to  strong  strain. 

2.3.  INITIAL  AND  BOUNDARY  CX)NDmONS 


(8) 


(9) 
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A  variable>daisity  temporal  vorticity  layer  of  finite  thickness  is  assumed  axtstO.  A  right- 
hamled  rectangular  coordinate  system  fx,y,z)  is  chosen  so  that  the  initial  vortici^  distribution  is 
aligned  with  the  positive  y-axis,  the  flow  being  uniftsm  in  the  streamwise  x-direction.  A  second- 
order  Gaussiatt  vtHticity  disttibotion  with  standard  deviation  cris  adopted  to  describe  the  variation 
of  vorticity  within  die  layer.  The  thickness  of  die  vorticity  layer,  (7,  and  the  fiee  stream  velocity  are 
ebt^n  as  lengdi  and  velocity  scales.  The  initial  vorticity  and  velocity  field  are  given,  respectively, 
by:  ca^xjO)  -  2l(a  =  (a^x,0)  =  0,  u(x,0)  =  erfizJa),  v(xjO)  =  Mxfi)  = 

0.  For  <7  =  i,  the  velocity  initial  flowfield  satisfies:  u(z  -^±oo}  s  ±  i.  Furthqniore,  die  layer  is 
assumed  periodic  in  die  streamwise  x-  and  qianwise  y-directions,  with  periodicity  lengths  Xx  and 
Xy  respectively,  and  is  unbounded  in  the  cross-stream  z-direction.  The  initial  density  distribution 
is  an  error-function  profile,  and  die  reference  (tensity  scale  is  chosen  so  that  the  low-density  fluid 
hasps  i. 

The  periodicity  boundary  conditions  introduce  some  difficulties  in  the  evaluation  of  the 
flowfield  and  of  its  gradient  since  we  must  consido*  the  image  system  of  the  transport  elements. 
This  system  yields  an  additional  term  which  must  be  added  to  the  velocity  induced  by  the  elements 
in  the  domain.  Unlike  the  2D  case,  this  term  may  not  be  deduced  from  a  potential  flow,  and  closed 
form  expressions  to  include  its  effect  are  not  known  (Ashurst  &  Meiburg,  1988).  In  the 
computations,  we  have  adopted  the  procedure  described  in  (Knio  &  Ghoniem,  1991),  which 
consists  of  confuting  directly  the  effea  of  the  dght  immediate  neighbors  of  die  elments,  and 
^iproximating  die  induced  velocity  of  the  images  which  lie  within  a  square  of  side  400  Xx  by 
imeipolation  on  a  fixed  grid.  White  previous  studies  have  only  included  the  cemtribution  of  a  small 
number  of  images  (Ashurst  &  Meiburg,  1988),  die  procedure  suggested  in  Knio  &  Ghoniem 
(1991)  is  prefored  because  it  yields  more  accurate  representations  of  the  velocity  and  velocity 
gradient  fields,  and  avoids  the  gemmuion  of  numerical  boundary  layers  at  the  spanwise  boundaries 
of  the  (tomain  (Inoue,  1989). 
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3.  RESULTS 


In  this  section,  we  focus  on  the  defonnation  of  the  primary  (essentially  two-dimensional) 
structure  due  to  the  3D  instabilities,  and  on  the  formation  of  secondary  structures  as  these 
instabilities  evolve  into  dieir  non-linear  stages.  We  attempt  to  unify  the  various  postulates  on  the 
origin  and  mechanisms  of  the  3D  motion.  By  comparing  our  results  to  available  experimental 
evidence  (Breidenthal,  1980, 1981;  Jimenez,  1983;  Jimenez  et  al.,  1985,  Bernal,  1981,  Lasheras 
&  Choi,  1988;  Lasheras  et  al.,  1986)  and  to  numerical  solutions  (Ashurst  &  Meiburg,  1988, 
Metcalfe  et  al.,  1987,  Lowery  et  al.,  1987;  Grinstdn  et  al.,  1989),  we  proceed  to  clarify  some  of 
the  aforementioned  issues,  and  point  to  strengths  and  deficiencies  of  previously  proposed  models. 
Another  objective  is  to  provide  a  base  solution  which,  in  turn,  is  used  to  highlight  the  baroclinic 
effects  in  the  variable-density  layer. 

3.1.  SYMMETRIC,  UNIFORM-DENSITY  SHEAR  LAYER 

A  temporal  shear  layer  with  streamwise  poiodicity  length  =  132,  which  matches  the 
wavelength  of  the  most  unstaUe  2D  mode  (Ghoniem  &  Krishnan,  1989)  and  spanwise  periodicity 
length  Ay  s  Acc/2,  which  lies  close  to  the  most  anq>lified  3D  nnode  of  the  translative  instability 
(Pienehumbert  &  Widnall,  1982),  was  computed.  The  initial  scalar  distribution  has  a  zero  mean 
and  a  unit  difference  across  the  layer.  The  shear  layer  is  initiaUy  discretized  among  elements 
distributed  on  a  grid  of  20x14x5  points  along  the  x-,y-,  and  z-directions  respectively.  Thus, 
conq>utational  elements  are  distributed  on  five  material  or  iso-scalar  surfaces.  The  selection  of  the 
number  of  material  surfaces  is  chosen  as  the  minimum  number  required  for  accurate  rqnesentation 
of  the  eigenfunctions  of  the  Kelvin-Helmholtz  instalnlity  (Ghoniem  et  al.,  1988).  The  core  radius 
of  the  smoothing  functions  is  chosen  so  tiiat  strong  overlap  among  the  cores  of  neighboring 
elements  is  ensured,  and  the  vorticity  of  tiie  elements  is  obtained  by  minimizing  the  integral  error 
between  assumed  and  discretized  vorticity  profiles  (Knio  &  Ghoniem,  1991).  In  the 
conq}utations,  we  set  ^  =  0^  and  the  time  step  At  -  0.1.  The  layer  is  initially  perturbed  in  both 
the  streamwise  and  spanwise  directions  by  disphudng  the  computational  elements  in  the  cross¬ 
stream  direction  using  sinewaves  of  amplitude  e  =  0.02Ax,  i.e.  by  u^g  the  transformation  Zi  z,- 
■f  esin(27CXilAx)  e  sin(2nyi/Ay)  (see  figure  1.) 

3.1.1.  DEFORMATION  OF  MATERIAL  SURFACE 

Figures  2  and  3  depict  perspective  views  of  the  iso-scalar  surfaces  initially  located  at  z  =  0, 
and  -1J2  respectively.  The  surface  initiaUy  lying  atz  =  0 represents  the  middle  surface  where 
most  of  the  vmticity  is  concentrated,  while  increasing  or  decreasing  the  value  of  z  corresponds  to 
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motion  towards  the  top  or  bottom  streams.  The  plots  are  generated  from  the  point  of  view  of  an 
observer  located  at  (48,24,48). 

Due  to  rollup  of  the  layer,  computational  elements  accumulate  in  the  core  which  forms 
during4.0  <  t<  8.0.  Previous  analysis  of  the  2D  solution  indicates  that,  for  the  present  amplitude 
of  the  streamwise  perturbation,  the  linear  stages  of  the  evolution  of  the  primary  2D  instability  end 
between  t  =  4.0  and  8.0,  followed  by  roll-up.  The  amplitude  of  the  spanwise  perturbation  remains 
small  for  t  <  8.0,  i.e.  its  an^lificatitHi  is  essentially  suppressed  during  the  growth  of  the  2D  mode. 

The  growth  of  the  eddy  core  in  the  mid-section  of  die  domain  continues  while  its  spanwise 
waviness  anqilifies.  The  amplitude  of  the  spanwise  perturbation  increases  significantly  along  the 
core,  an  indication  of  the  evolution  of  the  translative  instal^ty  (Kerrehumbert  &  Widnall,  1982). 
This  uneven  axial  displacement  of  the  spanwise  core  is  accompanied  by  an  out-of  phase 
deformation  of  the  braids  under  the  influence  of  the  streamwise  vorticity  generated  within  the 
cores.  The  growth  of  the  translative  instability  is  shown  by  plotting  the  row  of  elements  initially 
aligned  along  the  core  centerline  in  figure  3b.  Vortex  lines  aligned  with  the  axis  of  the  spamvise 
core  sufier  a  mild  net  deformation  in  the  streamwise  direction.  The  evolution  of  the  spanwise  core 
instability  occurs  such  that  vortex  lines  constantly  re-align  with  die  direction  of  strain  while  being 
stretched  along  their  axial  directicMi. 

The  fluid  motion  along  the  braids  is  illustrated  in  figure  3c.  The  stretching  of  the  braids 
leads  to  the  intensification  of  streamwise  vorticity  produced  as  the  braids  are  deformed  by  the 
growth  of  the  "translative"  instability  along  the  core,  and  strained  by  the  2D  flowfield.  Vortex 
rods,  which  extend  throughout  the  braids  and  are  wrapped  around  the  spanwise  cores,  form  as 
streamwise  vorticity  tolls  into  coherent  eddies.  The  material  surfaces  spin  around  the  streamwise 
axes  of  these  eddies,  which  are  located  at  the  streamwise  boundaries  and  middle  of  the  domain. 
This  motion  is  acconqianied  by  die  thinning  of  the  strip  in  the  region  separating  neighboring 
streamwise  vortices,  thus  producing  a  "hairpin"  vortex  configuration.  The  resulting  deformation 
of  the  flow  map  is  captured  by  the  mesh  refinement  algorithm  which  can  also  be  observed  in  the 
surface  plots.  Since  the  flow  is  inviscid,  the  division  and  change  of  shape  of  the  transport 
elements  describes  the  strain  field,  eqiecially  when  the  elements  approach  the  spanwise  core  or 
when  they  are  attracted  towards  die  axes  of  the  vortex  rods. 

The  adaptive  re^nse  of  the  numerical  scheme  is  illustrated  in  Table  1  where  the  surface 
area.  A,  and  of  the  number  of  elements,  N,  used  along  individual  material  surfaces  are  given. 
While  both  grow  rapidly  following  the  roll-up  of  die  vorticity  layer,  the  increase  in  the  number  of 
transport  elements  occurs  at  a  higher  rate  than  that  of  the  surface.  Though  the  middle  surface 
deforms  at  higher  rate  and  carries  a  larger  number  of  transport  elements  than  the  remaining 
surfaces,  small  zones  of  high  strain  exist  along  all  material  surfaces.  These  zones  necessitate  the 
introduction  of  new  elements,  an  effect  which  precedes  the  severe  deformatitm  of  the  surfaces. 
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Table  1 

Normalized  surface  area.  A,  and  number  of  elements,  N,  for  the  individual  material  layers. 


t-4.0 

t  =  8.0 

t=  12.0 

t=  16.0 

Material  Layer  1 
Location:  z- -132 

A  =  1.03 
N=1.00 

A  =  1.12 

N  =  120 

A  =  1.62 

N  =  1.68 

A  =326 
N  =  3.80 

Material  Layer.  2 
Location:  z  =  -0.66 

A  =  1.04 
N=1.00 

A  =  122 

N  =  133 

A  =231 

N  =  234 

A  =  4.79 
N  =  7.05 

Material  Layer  3 
Location:  z  =  0.0 

A  =  1.04 
N=1.00 

A  =  1.43 

N  =  1.49 

A  =  2.71 

N  =  324 

A  =5.41 
N  =  933 

Material  Layer  4 
Location:  z  =  0.66 

A  =  1.04 

N  =  1.00 

A  =  122 
133 

A  =231 

N  =  2J4 

A  =  4.79 
N  =  7.05 

Material  Layer  5 
Location:  z=  132 

A  =  1.03 
N=1.00 

II  II 

o 

A  =  1.62 
N=1.68 

A  =  326 
N  =  3.80 

Figure  4  shows  cross-sections  through  all  die  material  (con^utational)  surfaces  at  the  time 
the  computations  are  stopped,  t  s  1S.0.  (As  indicated  below,  the  structure  of  the  flow  field  does 
not  vary  appreciably  for  t  ^  16  due  to  saturation  of  the  instabilities.  This  additional  frame  is 
provided  to  emphasize  this  feature  of  the  confuted  flow;  the  subsequent  discussion  will  be  limited 
to  results  obtained  for  t  ^  16).  We  take  streamwise  cross-sections  through  the  braid  and  core  in  the 
planes  located  at  jc  s  2.0  (d)  andjc  =  6.6  (c),  respectively,  and  spanwise  sections  along  the  planes 
y  =  /.6  (b)  and y  ^33  (a).  Small  circles,  whose  radii  are  chosen  smaller  than  the  core  radius,  are 
drawn  to  mark  the  intersection  points  with  the  transport  elements.  The  streamwise  sections,  (c) 
and  (d),  show  how  the  streamwise  rods,  resulting  from  the  roll-up  of  the  braids,  give  rise  to  the 
formation  of  the  mushroom  structures  (Bernal  &  Roshko,  1986;  Lasheras,  Cho  &  Maxworthy, 
1986;  Lasheras  &  Choi,  1988).  They  also  depict  how  the  extension  of  these  rods  around  the 
spanwise  cores  results  in  the  establishment  of  a  double  mushroom  structure  (Lasheras  &  Choi, 
1988).  The  results  show  that  the  deformafion  of  the  material  surfaces  within  the  core  is  due  to  the 
combined  effect  of  three  rows  of  streamwise  vortex  structures:  two  rows  resulting  from  the 
extension  of  vcntex  rods  towards  the  core,  and  a  third  generated  by  the  deformation  of  the  core 
itself  under  the  action  of  the  translative  instability.  This  is  shown  schematically  in  figure  S. 

The  spanwise  sections  illustrate  the  effect  of  the  translative  instability  on  the  cross-stream 
position  of  the  core  and  the  shq)e  of  its  cross-section.  Figure  4  shows  that  the  core  is  pushed 
upwards  and  in  the  flow  direction  of  the  top  stream  in  the  "left”  half  of  the  domain,  0<y<  Xy/2, 
while  it  suffers  an  antisynunetric  deformation  in  the  other  half,  as  can  be  seen  in  figures  2  aiKl  3. 
The  core  loses  its  symmetry  at  most  spanwise  stations,  as  computational  elements  migrate  in  the 
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direction  opposite  to  that  of  the  core  translation.  The  distribution  of  material  particles  in  the  plane  y 
~  3.3,  which  intersects  the  central  streamwise  vortex  rod,  shows  that  the  braids  significantly 
thicken  at  tiiis  critical  spanwise  location  by  entraining  irrotational  fluid  from  one  side  of  die  layer  to 
the  other.  This  is  be  verified  by  simultaneously  examining  the  streamwise  cut  through  the  plane  x 
=  6.6,  which  illustrates  the  entrainment  of  die  mushrooms  around  the  spanwise  core. 

3.1.2.  VORTICITY  AND  SCALAR  DISTRffiUnONS 

The  motion  of  die  material  surfaces  follows  the  evolutim  of  the  vorticiQr  field.  This  motion 
establishes  entrainment  patterns  within  the  shear  layer.  Both  are  displayed  in^gures.  6-9  in  the 
form  of  vorticity  and  scalar  contours  plotted,  respectively,  on  two  ^anwise  sections,  y  =  3J,  and 
5.0,  and  on  two  streamwise  sections,  x  =  2.0  and  6.6.  The  vorticity  and  scalar  contours, 
generated  using  diffoent  techniques,  are  generated  at  times  t  =  4.0. 8.0, 12.0,  and  16.0.  Vorticity 
contours  are  generated  by  using  die  core  smoodiing  functior  <;  to  compute  the  vorticity  on  a  mesh  of 
40x40  points  and  processing  the  data  with  the  NCAR  contouring  software.  A  consideraUe  amount 
of  smoodiing  is  introduced  in  this  procedure,  and  the  resulting  plots  are  mainly  used  to  deduce  the 
large-scale  features  of  the  vorticity  field.  A  diffnent  approach,  m  which  shaded  areas  of  constant 
scalar  concentration  are  generated  by  interpolating  the  discrete  scalar  values  on  a  13Sxl35-cells 
grid,  is  adopted  in  the  representation  of  scalar  distribution. 

A.  Early  Stages  of  the  3D  Motion 

At  early  stages,  and  until  t  =  8.0,  the  vorticity  and  scalar  contours  plotted  at  various 
spanwise  sections  are  similar,  indicating  that  the  growth  of  the  3D  perturbations  is  suppressed 
during  the  early  stages  of  the  2D  instability  and  that  the  spanwise  vorticity  remains  essentially 
uniform  across  the  layer  (Metcalfe  et  aL,  1987).  Weak  streamwise  structures  which  change  their 
form  between  different  streamwise  stations  develop,  but  have  not  yet  gained  enough  strength  to 
alter  the  flow  significantly.  Mean^i^ule,  a  single  row  of  counter-rotating  vortices  is  found  to  repeat 
itself  at  all  streamwise  locations  of  the  domain.  These  structures  are  generated  by  tilting  of  the 
vortex  lines  which,  at  /  =  0,  do  not  possess  a  streamwise  vorticity  conqwnent  The  tilting  of  the 
vortex  lines  into  the  streamwise  ditection(s)  leads  to  die  creation  of  zones  of  alternating  streamwise 
vorticity  whose  locations  and  signs  follow  the  shape  of  sinewave  perturbation.  For  r  >  4,  the 
streamwise  vorticity  is  slowly  intensified  under  the  2D  strain  field,  producing  higher  values  in  the 
braids  of  the  eddy.  With  the  roll-up  of  the  vorticity  layer  and  the  formation  of  a  spanwise  core,  the 
edges  of  the  core  are  stretched  up  and  down  towards  the  free  streams,  giving  rise  the  top  and 
bottom  rows  of  counter-rotating  vortices  which  appear  in  figure  8a.  These  two  rows  are  separated 
by  a  third,  which  appears  as  a  small  circle  whose  vorticity  is  of  the  opposite  sign  to  the  previous 
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two.  As  shown  before,  the  latter  is  generated  as  a  result  of  the  growth  of  perturbations  on  the  core 
itself  by  the  translative  instat^ty  mechanism. 

B.  Late  Stages  of  3D  Motion 

For  8.0  <  t  <  12.0,  the  flowfield  suffers  a  rapid  transition  to  3D  motion  leading  to  an 
intensification  of  the  streamwise  vordcity  in  the  braids.  The  total  circulation  of  the  streamwise 
vortices,  =  fla>xl  d4  in  the  plane  jc  =  2.0,  increases  from  Fio^S)  =  2.617  to  Fto^n)  =  4.492. 
The  streamwise  vordcity  grows  under  the  acdon  of  the  strong  strain  exerted  by  the  large  spanwise 
cores  in  the  neighboriiood  of  the  stagnadon  "lines"  which  anchor  the  braids  (Lasheras  &  Choi, 
1988).  Meanwhile,  the  deformadon  of  the  core,  which  is  attributed  to  the  growth  of  die  transladve 
instability,  changes  the  alignment  of  the  vordcity  from  predominandy  spanwise  into  spanwise  and 
streamwise  components.  Although  the  maximum  value  of  streamwise  vordcity  occurs  within  the 
braids,  the  streamwise  eddies  generated  by  the  core  deformadon  have  higher  total  circuladon.  At  t 
=  12.0,  the  middle  row  of  streamwise  voidces  in  figure  8a  accounts  for  65%  of  the  total  circulaticm 
in  the  plane  of  the  core.  This  is  because  the  growth  and  maturadon  of  the  primary  2D  instability, 
which  precede  the  3D  modon,  force  the  migration  of  the  spanwise  vordcity  from  the  thinning 
braids  into  the  core. 

The  total  streamwise  circulation  in  the  plane  dividing  the  core,  x  s  6.6,  used  as  a  measure 
of  the  3D  effects  in  the  flow,  is  shown  in  figure  10.  It  confirms  the  early  observation  that  3D 
effects  are  small  during  the  linear  stages  of  the  2D  instability  and  grow  rapidly  after  roll-up.  The 
behavior  of  the  curve  changes  from  an  algetxaic  growth  for  r  <  9,  to  an  exponential  growth 
between  9  <t  <  13.  At  later  stages,  a  non-linear  regime  characterized  by  a  drop  in  the  rate  of 
increase  of  the  total  circuladon  is  observed.  This  is  expected  since  the  strain  field  induced  by  the 
spanwise  eddy  leads  to  continuous  intensification  of  the  streamwise  vordcity.  For  t  >  16.0,  no 
qualitative  changes  in  the  structure  of  the  vordci^  and  scalar  fields  is  observed,  an  indication  that 
the  instabilities  tend  to  saturate  (Metcalfe  et  al.,  1987). 

At  late  stages,  the  streamwise  vortices  induce  a  strong  secondary  motion.  As  previously 
indicated,  this  motion  can  be  easily  analyzed  in  the  braids  of  the  eddy  where  scalar  mushroom 
structures  are  generated.  However,  the  scalar  distribution  is  more  conq)lex  in  the  cote,  where  the 
flow  is  under  the  combiiied  influence  of  spanwise  and  streamwise  vortices.  The  top  and  bottom 
mushrooms,  which  form  due  to  the  roll-up  of  the  braid  vordcity  and  are  identified  in  figure  8b, 
originate  in  the  braids  and  are  then  entrained  towards  die  cores.  Near  the  axis  of  the  core,  the 
scalar  distribution  is  affected  by  the  flowfield  induced  the  three  rows  of  alternating  streamwise 
vortices  shown  schematk;ally  in  figure  5.  The  superposition  of  the  fields  of  tiiese  vortices,  whose 
axes  are  deformed  under  the  action  of  the  translative  instability,  leads  to  the  genoation  of  "W- 
shaped"  scalar  structures. 
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The  spanwise  vordcity  contours  at  r  =  12.0,  although  similar  to  those  encountered  in  a  2D 
flow,  exhibit  a  more  compact  core  than  in  2D  due  to  stretching  along  the  axis  of  the  core.  The 
stretching  of  spanwise  vorticity  is  accompanied  by  a  non-uniform  deformation  of  the  cores  at 
different  spanwise  locations.  The  corresponding  variation  of  the  spanwise  vorticity  is  depicted  in 
the  last  frames  of  figures  6  and  7.  The  core  is  shifted  downwards  and  in  the  direction  of  the 
bottom  stream  for  Xy/2  <y<Xy,  while  it  suffers  an  antisymmetric  displacement  for  0  <  y  <  Xy/2 
(see  also  figures  2-3).  This  confirms  the  results  of  the  linear  stability  theory  of  perturbed  vortex 
cores  which  predicts  a  "translative"  instability  of  the  cores  in  the  manner  described  above 
(Pierrehumbert  &  Widnall,  1982).  We  also  note  that  the  point  of  maximum  vorticity  within  the 
core,  which  moves  in  the  direction  opposite  to  diat  of  the  outer  boundaries  of  die  core,  no  longer 
coincides  with  its  geometric  center.  This  configuration  resendiles  that  observed  in  the  evolution  of 
the  eigenfunctions  of  the  linear  stability  problem  of  vortex  rings,  which  also  predicts  a  similar 
behavior  for  any  locally  curved  vortex  filament  (Widnall  and  Tsai,  1977).  We  recall  diat  the  most 
amplified  mode  of  the  translative  instability  mechanism  is  characterized  by  an  eigenfunction  which 
changes  sign  within  the  vorticity  core.  This  property  also  arises  in  the  linear  stability  of  vortex 
rings  as  a  necessary  condition  for  eigenfunction  instability.  In  both  cases,  the  amplification  of  the 
instability  forces  the  migration  of  the  "inner"  c<xe  in  the  direction  opposite  to  die  motion  of  its  outer 
boundaries. 

This  mechanism  ^ipears  to  be  connected  to  convective  currents  widiin  the  core  and  not  to  a 
uneven  vorticity  stretching.  This  is  verified  by  inflecting  die  scalar  distritnition  in  the  same  cross- 
sections  which  shows  that  the  scalar  field  follows  a  similar  redistribution  and  loses  its  symmetry. 
This  motion  leads  to  preferential  entrainment  of  irrotational  fluid  fivm  the  free  streams.  The 
section  of  the  core  which  is  displaced  upwards,  0<y<  Xy/2,  entrains  more  fluid  from  the  bottom 
stream,  while  the  section  which  is  pushed  downwards  consists  mainly  of  the  top-stream  fluid. 
This  form  of  preferential  entrainment  resembles  that  reported  experimentally  (Bernal  &  Roshko, 
1986). 

The  symmetry  of  die  vorticity  and  scalar  distributions  is  preserved  at  the  spanwise  mid¬ 
section  of  the  domain,  since,  as  predicted  by  die  themy,  the  curvature  of  the  core  vanishes  at  that 
location.  Nevertheless,  this  plane  is  of  interest  since  it  intersects  the  streamwise  vortex  rod 
centered  in  the  domain  (see  figures.  2-3).  This  rod  appears  in  the  form  of  a  "tongue"  of  negative 
vorticity  which  extends  duough  the  braids  to  the  top  and  bottom  edges  of  neighboring  c<»es  (figure 
6a,  r  =  16.0).  At  this  plane  of  zero  curvature,  the  action  of  the  translative  instability  is  manifested 
by  the  presence  of  two  vorticity  maxima.  A  similar  vorticity  distribution  is  obtained  in  unstable 
vortex  rings  at  azimuthal  stations  where  the  curvature  of  the  axis  of  the  core  vanishes  (Knio  & 
Ghoniem  (1988).  In  the  shear  layer,  these  stations  lie  within  the  planes  y  -0,  Xy/2,  and  Xy, 
while,  in  the  vortex  ring,  the  local  curvature  vanishes  whenever  the  curvature  induced  by  the 
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growth  of  azimuthal  bending  waves  cancels  that  of  the  undisturbed  vortex  ring.  The  resemblance 
between  the  two  cases  bears  important  implications  regarding  the  implementation  of  numerical 
schemes  and  the  naodeling  of  the  vorticity  layer.  The  numerical  study  in  Knio  &  Ghoniem  (1990) 
shows  that  the  dynamics  within  concentrated  vortices  may  not  be  properly  predicted  unless  a 
sufficiently  large  number  of  con^utational  elements  is  used  to  discretize  the  the  vortex  cores.  As  a 
result,  the  computations  coukl  have  missed  or  spuriously  predicted  the  evolution  of  the  translative 
instability,  had  we  chosen  to  simulate  tiie  vortktiQr  layer  by  distributing  tire  transput  elements  on  a 
single  material  layer,  or  to  model  the  layer  as  a  thin  vortex  sheet,  as  in  Ashurst  &  Meiburg  (1988) 
and  Inoue  (1989). 

3.1.3.  EbmiAINMENT  ENHANCEMENT 

The  development  of  the  3D  instabilities  promotes  shear  layer  entrainment  (Knio  & 
Ghoniem,  1991).  To  quantify  this  effect,  the  shear  layer  entrainment  is  measured  by  introducing 
an  "eddy  size"  parameter,  5,  in  2D  and  3D.  In  2D,  the  eddy  size  is  defined  by  the  measure  of  the 
region  enclosed  between  the  surfaces  where  the  normalized  scalar  first  deviates  by  3%  from  the 
corresponding  free  stream  value.  This  yields  a  local  height  of  the  eddy,  Z2d(x),  which,  when 
integrated  over  the  streamwise  length  of  the  domain,  pves  a  mean  eddy  size, 

Std  -  I  Z2d(x)  dx 

JO  (10) 

According  to  this  definition,  the  eddy  consists  of  the  union  of  the  rotational  fluid  and  the 
irrotational  fluid  trapped  between  the  braids  and  the  core  (see  figure  11).  In  3D,  two  similar 
definitions  are  considered.  The  first  is  obtained  1^  measuring  the  height  of  the  eddy  at  each 
spanwise  and  streamwise  location,  Zjofx.y),  and  then  integrating  over  the  spanwise  and 
streamwise  periodicity  length  to  get  an  eddy  size, 

JrA,  fXy 

I  I  ZiDix^y)  dy  dx 

0  Jo  (11) 

While  this  definition  is  a  natural  extension  of  diat  in  2D,  the  contribution  of  the  irrotational  fluid 
trapped  by  the  mushroom  structures  is  neglected  in  the  averaging  process.  To  account  for  this 
additional  entrainment  mechanism,  we  define  a  third  measure,  S'jd,  by: 
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(12) 


53D  =  I  23x)(x)  dx 
JO 

where  Z'si^x)  -  maxy(Z3D(x,y)). 

The  eddy  size  is  normalized  by  the  spanwise  and  streamwise  periodicity  lengths  in  3D 
computations  and  the  streamwise  periodicity  length  in  2D,  so  diat  die  resulting  values  represent 
an  average  thickness  of  the  scalar  distribution.  The  entrainment  enhanconent  is  shown  in  figure  12 
by  comparing  the  eddy  size  for  2D  and  3D  computations.  Widi  the  roll-up  of  the  layer,  t  ~  8.0, 
entrainment  curves  start  to  grow  with  the  curves  in  3D  acquiring  a  hi^er  growth  rate.  The 
deviation  in  the  behavior  of  the  2D  and  3D  solutions  coincides  widi  die  3D  transition  depicted  in 
figure  10.  Thus,  the  transition  to  3D  motion  is  accompanied  by  an  "oitrainment  transition."  By 
the  end  of  the  simulation,  the  total  entrainment  increases  by  75%  over  its  2D  counterpart, 
while  comparison  of  S'3d  tutd  indicates  that  the  formation  of  the  mushroom  structures 
contributes  significandy  to  entrainment  enhancement 

The  eddy  size  parameters  are  also  used  to  quantify  the  preferential  entrainment  of  fluid  at 
various  sections  of  the  domain.  This  form  of  entrainment  reverses  itself  every  half  spanwise 
wavelength  so  that  die  composition  of  the  eddy  does  not  favor  eidier  firee  streartt  However,  in  the 
region  0  <  y  <  Ay2,  preferential  entrainment  of  lower  stream  fluid  is  observed.  Preferential 
entrainment  may  be  estimated  by  subdividing  die  integrals  in  Eqs.  (10-12)  over  the  regions  defined 
by:  s  >  Jav.  and  r  <  Jav,  thus  yielding:  (830)*.  (830)’,  (8'3d)*^  and  (5'jp)'.  At  t  =  16.0,  the 
entrainment  ratios  in  the  area  0  <  y  <  kyl2.  =  1.787,  and  (8’3d)-I(8'3d)^  =  3.487. 

These  ratios  are  reversed  in  the  area  Xyf2  <y  <  Xy,  so  diat  unit  net  entrainment  ratios  are  obtained. 
The  difference  between  the  two  ratios  is  a  manifestation  of  the  role  of  streamwise  vortices  in 
inducing  preferential  entrainment  patterns  (Bernal  &  Roshko,  1986). 

3.1.4  DISCUSSION 

The  results  of  the  conputations  show  that  the  evolutitm  of  die  shear  layer  from  a  perturbed 
steady  state  using  monochromatic  3D  distuibances  consists  of  diree  stages:  (1)  In  the  first  stage, 
an  essentially  2D  growth  of  the  perturbations,  in  die  form  of  linear  anplification  of  the  Kelvin- 
Helihholtz  instability  modes,  is  observed.  During  this  phase,  all  3D  activity  is  suppressed.  This  is 
followed  a  non-linear  2D  growth  of  die  fundamental  mode,  leading  to  toU-ip  and  die  formation 
of  a  concentrated  panwise  eddy  core.  (2)  The  roU-ip  is  acconpanied  a  rapid  growth  of  the  3D 
modes,  in  the  form  of  a  deformation  of  the  panwise  eddy  and  an  instability  in  the  braids.  (3)  The 
3D  modes  undergo  a  non-linear  growth  which  results  in  the  redistribution  of  die  streamwise 
vorticity  into  vortex  rods  and  in  die  generation  of  the  scalar  mushroom  structures. 
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While  the  2D  instability  is  relatively  simple  to  describe,  the  3D  motion  in  the  cores  and 
braids  of  the  vortex  structures  are  more  involved  The  stability  of  the  rolled  layer  was  examined  by 
Pierrehumbert  and  Widnall  (1982)  using  a  periodic  array  of  Stuart  vortices.  They  found  that  this 
array  was  linearly  unstable  to  a  mode  they  called  the  translative  instability,  which  leads  to  a 
deformation  of  the  spanwise  eddies  according  to  the  strain  field  induced  by  the  eddy  and  by  its 
image  vortices.  The  instability  appearing  in  the  braids  was  studied  in  the  work  of  Corcos  &  Lin 
(1984),  who  proved  the  fundamental  nature  of  the  instability  of  a  strained  stteamwise  vorticity  and 
showed  that  the  growth  of  this  mode  is  suppressed  during  the  anq>lification  of  the  2D  instability. 
The  translative  instability  is  fundamentally  different  from  die  Corcos  mechanism  since  the  braids  of 
the  spanwise  vortices,  observed  in  the  case  of  a  shear  layer,  are  not  well  represented  by  the  Stuart 
vortices.  Thus,  the  two  instabilities  differ  in  fomL  The  translative  mode  is  an  instability  of  the 
spanwise  cores,  that  is  of  large  concentrated  vortices  subject  to  strain  normal  to  their  axes,  while 
the  Corcos  mechanism  predicts  the  instability  of  streamwise  vorticity  when  subjected  to  strong 
extensional  strain. 

The  absence  of  toaids  does  not  preclude  the  formation  of  vortex  rods  by  the  strain  field 
which  may  cause  the  vorticity  widiin  the  core  to  migrate  outwards  preferentially  (Grinstein  et  aL, 
1989).  On  the  other  hand,  the  ability  of  the  streamwise  vortex  rods  generated  by  the  Corcos 
mechanism  to  impart  a  core  deformation  similar  in  shqie  to  that  obtained  by  the  development  of  the 
translative  instability  (Corcos  and  Lin,  1984)  complicates  the  task  of  separating  the  role  of  these 
two  instabilities.  This  difficulty  has  led  researchers  (Ashurst  &  Meiburg,  1988;  Lasheras  &  Choi, 
1988)  to  en^hasize  the  importance  of  a  non-linear  interaction  between  the  ^>anwise  cores  and  the 
streamwise  vortices.  Our  results,  however,  indicate  that  the  translative  instability  within  the  cores 
grows  simultaneously  with  the  intensification  of  the  rods.  Sirrce,  as  previously  mentioned,  a  mild 
straining  of  the  streamwise  vorticity  is  observed  prior  to  the  roll-up  of  the  spanwise  vorticity,  and 
since  the  translative  instal^ty  is  a  linear  instability  mechanism,  the  importance  of  such  an 
interaction  should  be  dowrqrlayed. 

The  formation  of  the  streamwise  vortex  rods  leaves  its  mark  on  the  flow  in  the  form  of 
scalar  mushroom  structures.  The  translative  instability,  on  the  other  hand,  is  rrumifested  by  an 
asymmetric  q)anwise  vorticity  distribution  within  the  core,  which  resembles  that  predicted  by  the 
growth  of  unstable  eigoifunctions  of  the  Widnall  instability  of  vortex  rings  (Widnall  &  Tsai, 
1977).  The  alternating  preferential  entrainment  patterns  which  result  firom  the  amplification  of  the 
instability  mode  compound  the  difficulty  in  experimentally  visualizing  the  scalar  structures 
mentioned  by  Corcos  &  Lin  (1984),  since  uniform-density  spatially-developing  shear  layers 
entrain  more  fluid  fi'om  the  faster  stream  (Dimotakis,  1986).  Despite  these  difficulties,  a 
deformadon  consistent  with  the  develc^nnent  of  the  translative  instability  was  inferred  by  Jimenez 
(1983)  and  by  Bemal  &  Roshko  (1986)  based  on  spanwise  correlations  of  velocity  fluctuations. 
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The  analysis  of  the  evolving  vorticity  structures  using  numerical  simulation  makes  a 
valuable  contribution  since  direct  and  detailed  measurements  of  the  vorticity  Held  are  usually  not 
possible  or  extremely  cumbersome.  Such  difficulties  are  primarily  encountered  in  experimental 
studies  in  which  velocity  correlations  (Browand  &  Troutt,  1980;  Jimenez,  1983;  Wygnansky  et 
al.,  1979),  passive  scalar  or  dye  techniques  (Bernal  &  Roshko,  1986;  Jimenez  et  al.,  1985),  or 
low  heat  release  chemical  reactions  (Breidenthal  1981;  Lasheras  et  al.,  1986;  Lasheras  &  Choi, 
1988)  have  been  used  as  substitute  tools  for  deducing  the  topology  of  the  vorticity  field.  Our 
results  suggest  that  such  efforts  should  be  conducted  with  great  care,  because  a  large  number  of 
mechanisms  contributes  to  their  formation.  The  conq)uted  results  reveal  some  of  the  difficulties 
arising  from  the  lack  of  an  accurate  knowledge  of  the  vorticity  field.  In  particular,  it  is  shown  that 
the  deformation  of  material  surfaces,  visualized  fw  instance  by  injecting  marker  particles  in  one  of 
the  fluid  streams,  may  not  be  sufficient  to  fully  determine  the  corresponding  vortical  structures, 
especially  if  one  cannot  a  priori  locate  a  given  material  surface  with  respect  to  the  layer  of  highest 
spanwise  vorticity.  The  rapid  variation  of  the  strain  field  around  the  latter  causes  substantially 
different  deformations  of  adjacent  layers,  thus  preventing  immediate  correlations  with  the 
underlying  vortical  structures.  On  the  othm  hand,  the  visualization  of  v<^city  structures  by  means 
of  the  products  of  a  unity-stoichiometry,  low-heat-telease  chemical  reaction  may  remove  some  of 
these  difficulties,  since,  as  observed  by  Knio  and  Ghoniem  (1992),  the  products  of  reaction  are 
always  entrained  into  zones  of  hi^  vorticity. 
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3.2  VARIABLE-DENSITY  SHEAR  LAYER 

Density  variation,  characterized  by  a  density  ratio  of  the  free  streams,  plays  an  important 
role  in  the  evolution  of  heterogeneous  shear  layers  created  by  the  merging  two  streams  of  unequal 
density  and  velocity.  It  has  been  observed  that  a  non-unity  density  ratio  alters  die  spatial  growth  of 
the  layer  and  influences  the  entrainment  induced  by  the  vortical  structures  embedded  therein,  even 
when  gravity  effects  are  weak  (Brown  &  Roshko,  1974;  Ho  &  Huene,  1984;  Dimotalds,  1989). 
In  chemically-reacting  flows,  density  variation  is  generated  by  heat  release  which  leads  to  the 
generation  of  zones  of  high  tenqierature  and  low  density.  Here  too,  the  presence  of  two  or  more 
zones  of  different  density  is  found  to  affect  the  stability  and  development  of  the  flow.  This  effect 
depends  on  both  the  details  of  the  density  and  vorticity  distributions,  as  it  may  constitute  a 
stabilizing  or  a  destabilizing  mechanism  (Riley  &  McMurtiy,  1989;  McMurtry  et  al.  [47];  Ghoniem 
&  Krishnan,  1988). 

In  this  section,  the  effect  of  weak  density  variation  on  3D  instability  is  investigated  by 
computing  the  onset  of  3D  motion  in  an  inconqiressible,  variable-density  shear  layer.  While  this 
study  is  ultimately  motivated  by  the  desire  to  predict  high-heat-release  reactions  in  vortical  flows, 
the  simplified  model  allows  us  to  focus  on  the  dynamic  effects  of  baroclinic  vorticity.  In  a  3D 
flow,  this  stepwise  t4)proach  is  important  because  of  the  presem^  of  a  vorticity  stretching  term. 
One  case  with  small  density  ratio,  (hnaxIPmin  -  2,  is  considered.  To  avoid  a  re-initialization  of  the 
scalar  (density)  gradient  field,  we  use  the  same  discrete  scalar  gradient  values  of  the  previous  case 
but  set  pav  s  IJ.  Therefore,  in  the  high-density  tc^  finee  stream,  p^u  =  2,  and  in  the  low-density 
bottom  free  stream,  pnon  -  f  •  The  same  initial  perturbation  used  in  the  previous  case  is  applied, 
and  the  computations  are  pcrfcxmed  to  observe  the  growth  of  both  the  2D  and  3D  instabilities. 

The  linear  theory  of  variable-density  shear  layers  (Ghoniem  &  Krishnan,  1989;  Krishnan, 
1989)  shows  tiiat  the  wavelength  of  die  most  amplified  mode  of  die  2D  instability  deprads  weakly 
on,  and  the  growth  rate  are  almost  indqiendent  of  the  density  ratio,  while  its  phase  velocity  is 
strongly  varied  strongly  with  the  density  ratio.  Unlike  the  uniform-density  case,  the  most 
amplified  mode  in  the  variable-density  layer  evtdves  as  a  traveling  wave  moving  in  the  directirai  of 
the  high-density  stream  with  phase  speed  increasing  widi  the  density  ratio.  This  effect  has  been 
used  to  explain  the  difference  in  growth  rates  in  variable-density,  spatially-evolving  layers 
(Dimotalds,  1989, 1986;  Brown  &  Roshko,  1974,  Ghoniem  &  Krishnan  1989).  For  the  density 
and  vorticity  {vofiles  used  in  the  simulation,  and  a  density  ratio  of  2,  the  linear  stability  themy 
predicts  a  phase  speed  c  s  0.17.  We  note  that  the  stability  analysis  performed  by  Ghoniem  & 
Krishnan  (1989)  differs  slightly  from  the  eariier  investigation  of  Maslowe  Sc.  Kelly  (1971),  who 
showed  that,  in  the  limit  of  vanishing  Richardson  number,  density  stratification  tends  to  stabilize  a 
temporal  shear  layer.  The  differing  conclusions  are  due  to  the  assumed  doisity  profiles.  The 
analysis  of  Maslowe  and  Kelly,  which  was  primarily  concerned  with  aunospheric  and  oceanic 
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flows,  assumed  an  exponendally-vaiying  density  distribution.  This  distribution  is  contrasted  with 
the  error-function  profile  which  is  more  representative  of  high-speed  shear  flows,  where  small 
Richardson  numbers  are  more  frequently  encountered. 

The  effect  of  the  density  variation  on  the  linear  stability  of  3D  pertutbatitms  is  not  known  at 
present  3D  instability  of  variable-density  shear  layers  is  not  a  simple  extension  of  that  of  the 
uniform-density  case  since  one  must  deal  with  the  ackled  difficulty  of  formulating  a  steady  initial 
condition  of  the  stability  problem.  In  the  uniform-density  case,  Stuart  vortices  were  used  to 
approximate  the  flow  of  the  rolled  layer  (Pierrehumbert  &  Widnall,  1982).  The  existence  of 
similar  solutions  in  the  variaUe  density  case  is  coRq)licated  tty  the  convective  motion  of  the  vortices 
due  to  die  baroclinic  generation  of  vordcity,  as  observed  in  2D  simulations  (Ghoniem  &  Krishnan, 
1989).  As  a  result,  and  because  we  have  approached  the  study  of  3D  tenqioral  layers  with  the 
intention  of  investigating  deviations  from  the  2D  case,  we  were  content  with  keeping  the  same 
spanwise  periodicity  length  used  in  the  uniform-density  case.  This  pronqited  us  to  select  a  small 
value  of  die  density  ratio,  so  that  conqiarisons  with  the  preceding  results  are  justified. 

3.2.1.  DEFORMATION  OF  MATERIAL  SURFACES 

Figure  13  shows  the  evolution  of  a  material  surface  initially  lying  in  the  planes  z^O.  At 
early  stages,  the  motion  of  the  middle  layer  bears  strong  resemblance  to  its  counterpart  in  the 
uniform-density  case  (figure  2).  This  simUarity  is  expected,  and  in  agreement  with  the  results  of 
the  linear  stability  theory  which  predicts  almost  identical  growdi  rates  of  die  2D  perturbation.  The 
first  manifestation  of  the  convective  motion  of  the  instability  wave  is  observed  once  the  roU-up  of 
the  layer  is  completed  and  a  well-defined  eddy  core  is  formed  (/  >  8).  During  this  period,  the  3D 
perturbation  is  anqilified  causing  the  core  to  defmm  in  the  way  similar  to  that  predicted  by  the 
translative  instability  of  uniform-density  cores.  However,  the  (deformed)  axis  of  the  core  no 
longer  coincides  with  the  streamwise  mid-section  of  the  domain,  but  shifts  in  the  positive 
streamwise  x-  direction,  i.e.  in  the  direction  of  the  high-density  top  stream.  The  braids  suffer  a 
deformation  whose  shiqie  is  of  the  same  type  obsoved  in  the  uniform-density  case. 

The  total  surface  area  of  the  material  surfaces  and  the  number  of  tran^rt  elements,  shown 
in  Table  2,  exhibit  the  same  trends  described  in  the  uniform-density  case.  However,  the  motions 
of  surfaces  lying  on  the  top  and  bottom  sides  of  the  middle  layer  are  no  longer  similar  because  of 
the  asymnnetry  of  die  flowfield.  Prior  to  die  nutfuration  of  die  3D  instability,  t  <  12.0,  the  material 
surfaces  lying  on  die  low-doisity  side  are  deformed  at  higher  rate  than  those  located  in  the  heavier 
fluid  side  of  the  layer.  At  later  stages,  this  trend  is  reversed.  This  asymmetric  straining  of  the 
surfaces  is  related  to  the  asymmetry  of  the  spanwise  vortices  whose  rotation  resembles  the 
eccentric  ginning  of  egg-shaped  cores. 
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While  the  uneven  deformation  of  the  braids  is  readily  verified  at  later  stages,  the 
asymmetric  stretching  of  stieamwise  vordcity  affects  the  entrainment  patterns  in  die  neighborhood 
of  the  spanwise  cores  (Brown  &  Roshko,  1974;  Dimotalds,  1986,  1989;  Ghoniem  &  Krishnan, 
1989).  However,  the  motion  of  the  Lagrangian  mesh  cannot  be  used  to  deduce  this  effect  because 
the  vorticity  lines  in  a  vaiiable-denrity  flow  cannot  be  directly  identified  with  material  lines,  and  the 
defonnation  of  material  surfaces  is  not  sufficient  to  con^letely  determine  the  fate  of  the  vorticity. 
In  fact,  detailed  examination  of  the  motion  of  die  matoial  surfaces  only  indicates  that,  in  die  braids, 
streamwise  vorticity  re-organizes  into  rods,  leading  to  the  formation  of  the  mushroom  structures 
that  cover  the  eddy  cores,  and  that  die  spanwise  cores  propagate  in  an  asymmetric  flowfield. 

Table  2 

Normalized  surface  area.  A,  and  number  of  elements,  N,  for  the  individual  material  layers. 


t  =  4.0 

II 

Oo 

t=  12.0 

t=16.0 

Material  Layer  1 
Location:  z  =  -132 

A  =  1J03 
N^l.00 

'"-1 

II  II 

A  =  1.70 

N  =  2.00 

A  =  3.47 
N  =  4.75 

Material  Layer  2 
Location:  z  s  ^.66 

A  ^1.04 
N^l.00 

A  =  122 

N  =  1.40 

A  =  235 

N  =  2.84 

A  =  4.74 
N  =  7.75 

Material  Layer  3 
Location:  z  s  0.0 

A  ^1.04 
N^IJOO 

A  =  1.43 
N=135 

A  =  2.71 

N  =  3.65 

ii  II 

Material  Layer  4 
Location:  z  =  0.66 

A^1j04 

N^l.00 

A  =  121 
N=1.41 

A  =  222 

N  =  2.93 

A  =  4.78 
N  =  8.63 

Material  Layer  5 
Location:  z  =  132 

11  11 

A  =  1.12 

N  =  128 

A  =  137 

N  =  1.71 

A  =  3.11 
N  =  3.80 

Signs  of  die  influence  of  the  baroclinic  vortici^  and  the  associated  asymmetry  of  die  strain 
field  on  the  evolution  of  the  3D  paturbadon  appear  in  the  cross-sections  dirough  the  material 
surfaces.  These  sections,  shown  at  r  s  16.0  in  figure  14,  are  generated  on  a  streamwise  section 
translated  in  die  direction  of  the  heavy  stream  tox  s  93  in  order  account  for  the  convective  motimi 
of  the  eddy.  The  mean  streamwise  location  of  the  cme  is  first  estimated  by  translating  the  plane 
X1I2  ~  6.6  by  ct,  c  bdng  the  (diase  speed  of  the  2D  Kelvin-Helmholtz  wave  obtained  from  the 
linear  theory.  This  estimate  is  thoi  refined  by  considering  nei^boring  planes  on  both  sides  of  die 
planex  s  xj/2  +  ct.  It  is  found  diat  die  mean  streamwise  location  of  the  eddy  shifts  fiomx  =  6.6  at 
f  =  0,  to  X  =  73,  8.0,  8.6,  and  9.3  at  t  ^  4.0,  8.0,  12.0,  and  16.0  respectively.  Thus,  the 
convective  speed  of  the  eddy  is  closely  qiproximated  by  c  evoi  in  the  non-linear  stages  (Ghoniem 
&  Krishnan,  1989).  For  the  remaining  sections,  we  still  use  the  planes  y  ^  33  and  1.6  to 
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visualize  Ac  variations  along  the  span  of  the  layo',  and  the  plane  x  =  2.0  to  obtain  a  representative 
streamwise  section  of  the  braids. 

The  large-scale  features  of  the  instability  in  the  variable-density  flow  can  be  approximated 
by  those  found  in  die  previous  case  simply  by  accounting  for  the  convective  motion  of  the  eddy. 
However,  cross-sections  through  the  core  reveal  that  the  mushrooms  entrained  on  top  and  bottom 
of  the  core  are  not  similar.  The  top  mushroom,  lying  on  the  side  of  the  high-density  fluid,  is 
larger,  more  rounded,  and  less  developed  than  its  counterpart  on  the  bottom  side  of  the  eddy.  The 
concentration  of  computational  elements,  visualized  by  darker  areas  on  the  plots,  is  higher  in  the 
lower  mushroom  especially  near  the  axes  of  the  streamwise  vortex  rods.  Thus,  streamwise 
vorticity  is  higher  for  the  bottom  vortex  rods,  leading  us  to  expect  higher  rates  of  spinning  around 
their  axes,  and  a  significant  dqiaiture  from  the  entrainment  patterns  observed  in  the  previous  case. 

3.2.2  VORTICrrY  AND  DENSITY  FIELDS 

The  spanwise  structure  of  the  layer  is  shown  in  terms  of  the  spanwise  vorticity  and  density 
contours,  plotted  in  figures  IS  and  16  in  the  spanwise  planes  y  =  33  and  5.0,  respectively.  The 
streamwise  structure  of  the  layer  is  examined  by  considering  cross-sections  through  the  core  and 
the  braid.  The  core  cross-sections,  figure  17,  are  generated  in  the  y-z  planes  coinciding  with  the 
mean  streamwise  location  of  the  eddy.  We  use  the  fixed  streamwise  plane  x  =  2.0  to  cut  through 
the  braids,  the  results  being  shown  in  figure  18.  Results  are  shown  at  r  =  12.0  and  16.0,  i.e. 
following  die  growth  of  the  3D  modes. 

A  qualitative  smtilarity  between  the  uniform-  and  variable-doisity  flow  in  terms  of  the  Q^ies 
and  shapes  of  vortical  structures  that  are  formed  as  a  result  of  the  evolution  of  the  various 
instabilities  is  noticeable.  The  variable-density  layer  can  be  characterized  by  the  same  stages  of 
evolution  as  the  uniform-density  layer,  (i)  an  early  growth  of  the  Kelvin-Helmholtz  instability 
during  which  the  layer  remains  essentially  2D;  (ii)  a  non-linear  evolution  accompanied  by  the 
formation  of  a  spanwise  core  as  a  coherent  eddy,  the  onset  of  the  3D  undulation  along  its  axis,  and 
the  generation  of  streamwise  vorticity;  and  (iii)  a  maturation  of  the  translative  instability,  the 
redistribution  of  die  streamwise  vorticity  into  vortex  rods,  and  the  formation  of  scalar  mushroom 
structures.  The  differences  between  the  two  cases,  which  arise  due  to  the  baroclinic  generation  of 
vorticity,  are  summarized  in  die  following: 

(1)  The  evolution  of  the  Kelvin-Helmholtz  mode  is  modified  by  a  finite  phase  speed  of  the  waves 
in  the  direction  of  the  high-density  stream.  The  motion  of  die  waves  is  uniform  in  all  spanwise 
stations,  and  persists  into  (he  non-linear  stages. 

(2)  The  vorticity  field  loses  its  symmetry  as  a  result  of  the  vorticity  generated  by  the  baroclinic 
torque.  The  loss  of  symmetry  is  not  restricted  to  any  particular  spanwise  plane,  and  is  not  due 
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to  the  amplification  of  3D  modes  as  it  is  also  observed  in  2D  simulations  (Ghoniem  & 
Kiishnan,  1989). 

(3)  The  loss  of  symmetry  afiects  the  stieamwise  vorticity  via  its  asymmetric  strain  field.  For  t  < 
12.0,  the  layer  of  streamwise  vortices,  lying  on  the  side  of  the  high-density  stream,  is 
consideraUy  weaker  than  its  counterpart  on  the  bottom  side,  while  the  trend  is  reversed  at  later 
stages.  At  f  =  12.0,  Fto/^rbo,  =  0J36,  and  increases  to  ftopirbot  =  1.12  at  r  =  16.0.  Thus, 
the  strength  of  the  stieamwise  vortices  changes  on  both  sides  of  the  eddy  in  accordance  with 
the  deformation  of  material  surfaces.  This  implies  that,  for  12.0  <t<  16.0,  the  asymmetric 
straining  of  streamwise  vorticity  leads  to  the  establishment  of  spanwise  entrainment  patterns 
which  are  biased  towards  the  high-density  side.  This  bias  is  due  to  the  difference  in  strength  of 
the  streamwise  rods  wriyiped  on  opposite  sides  of  the  spanwise  core,  and  opposes  the  efiect  of 
the  streamwise  (2D)  entrainment  currents  v^ch  favor  low-density  fluid. 

(4)  The  loss  of  symmetry  is  aocompanied  by  a  net  asymmetric  entraiiunent  of  the  low-dmsity  fluid. 
We  distinguish  between  the  preferential  entrainment  associated  with  the  growth  of  3D  modes, 
and  the  entrainment  asymmetry  due  to  the  baroclinic  generation  of  vorticity.  The  asymmetric 
entrainment  of  the  low-density  fluid  combines  with  the  preferential  entraiiunent  of  low-density 
fluid  in  0  <  y  <  2^2,  produces  an  asymmetric  distribution  at  the  spanwise  mid-section  of  the 
domain,  and  counteracts  the  effect  of  the  3D  instabilities  in  Ay2  <y  <Xy.  At  r  s  16.0,  the 
enurainment  ratios,  (SiJo)’l(S$o)*  =  2541,  and  (S'jpMS'io)'*"  =  4.458  in  the  area  0  <y  < 
Xy/2,  while  the  net  entrainment  ratios,  (SsD)’f(S3D)*  -  1245,  and  (5'5£))'/(S'jr£))+  =  1.140. 
Thus,  the  spanwise  entrainment  patterns  induced  by  the  stieamwise  vortex  rods  lead  to  a 
reduction  of  die  asymmetry  in  the  entraiiunent  ratio. 

3.2.3  BAROCLINIC  VORTICrrY 

In  order  to  isolate  the  effects  of  the  density  variation  from  those  associated  with  vortex 
stretching,  the  distribution  of  the  baroclinic  torque  is  used.  Figure  19  shows  the  spanwise 
conqionent  of  the  baroclinic  torque,  in  the  planes  y  ^  33  and  5.0,  while  the  streamwise 
component  is  shown  in  figure  20  for  the  core  tuid  braid  sections.  At  early  stages,  t  <  8.0,  the 
spanwise  coiiqxinent  of  the  baroclinic  torque  is  coiuxntrated  in  two  zones  of  opposite  signs.  The 
initial  vorticity  of  the  layer  is  dqileted  in  the  left-hand  side  of  the  domain  and  enhanced  in  the 
remaining  part  Thus,  baroclinic  vorticity  inqiarts  an  asymmetry  on  the  vorticity  distribution  such 
that  the  pan  of  the  layer  di^laced  towards  the  high-density  fluid  is  weakened,  while  diat  pushed  in 
the  direction  of  the  low-density  stream  is  intensified.  As  suggested  in  Ghoniem  &  Krishnan 
(1989),  this  asymmetry  may  be  used  to  explain  the  origin  of  the  motion  of  the  Kelvin-Helmholtz 
mode. 
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Reduction  of  spanwise  vorticity  on  the  high-density  side  of  the  layer,  and  its  enhancement 
on  the  low-density  side  is  responsible  for  the  observed  difference  in  the  strength  of  the  streamwise 
vortices  between  the  top  and  bottom  rows  of  the  core.  As  shown  in  figure  20,  the  streamwise 
component  of  the  baroclinic  term  is  extremely  weak  during  the  linear  stages  of  the  primary 
instability.  Therefore,  the  asynuneny  in  the  streamwise  vorticity  distribution  must  be  due  to  the 
uneven  tilting  and  stretching  of  layers  of  varying  strengths  so  that  density  variation  effects  are 
primarily  felt  dirough  changes  in  the  spanwise  vorticity. 

At  later  stages,  however,  baroclinic  torques  becomes  strong  enough  to  directly  affect  the 
evolution  of  the  streamwise  vorticity.  Density  variation  does  not  lead  to  a  net  intensification  or 
weakening  of  the  streamwise  vortices  in  the  braids  since  the  baroclinic  term  changes  sign  within 
each  streamwise  eckly.  Nevettiieless,  the  baroclinic  tom  is  distributed  in  such  a  way  as  to  weaken 
the  top  parts  of  the  streamwise  eddies  and  to  strengthen  their  bottom  parts,  leading  to  a  downward 
drift  of  the  vortex  rods.  On  the  other  hand,  in  the  core,  baroclinic  torques  contribute  to  the 
asymmetry  between  the  top  and  bottom  rows  of  streamvidse  vortices.  While  the  middle  row  of 
vortices  is  affected  in  a  similar  way  as  that  observed  in  the  braids,  figure  20  indicates  that 
baroclinic  torques  tend  to  strengths  the  vorticity  of  die  botttmi  row  at  the  expense  of  the  top  row. 

Thus,  density  variation  plays  two  different  roles  in  the  development  of  3D  form  of  the 
layer.  At  the  early  stages,  it  generates  an  an  asymmetric  strain  field  by  inq)arting  a  convective 
motion  on  the  core.  At  the  later  stages,  it  redistributes  the  vorticity  within  the  core.  The 
development  of  the  instability  modes  in  die  variable-density  shear  layer  highlights  the  inqxnrtance 
of  tile  asymmetry  of  the  flow  and  strain  fields.  In  the  interpretation  of  the  origin  of  the  coriqilicated 
structures  associated  with  the  3D  effects,  the  influence  of  strain  and  density  gradient,  which  are 
respectively  taken  into  account  in  the  equation  of  motion  through  the  vorticity  stretching  and 
baroclinic  production  terms,  may  be  hard  to  distinguish  in  the  results.  To  facilitate  this  task,  we 
consider  the  case  of  a  uniform-density  shear  layn  with  an  asynunetric  vorticity  distribution  at  t  =  0, 
and  posqxMie  further  discussion  until  results  of  the  latter  case  are  analyzed. 
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3.3  UNIFORM-DENSITY.  ASYMMETRIC  SHEAR  LAYER 

Shear  layers  with  asymmetric  vorticity  profiles  are  frequently  encountered.  Typically, 
these  layers  are  formed  following  the  merging,  downstream  of  the  splitter  plates,  of  boundary 
layers  of  unequal  thicknesses  and  opposite  signs  of  voiticity.  The  velocity  profile  associated  with 
the  asymmetric  vorticity  distribution  thus  produced  can  be  modeled  as  the  superposition  of  a 
symmetric  velocity  profile  induced  by  an  idealized  symmetric  vorticity  distribution,  and  a  wake 
component  biased  towards  the  low-velocity  stream.  As  shown  by  Koochesfahani  &  Fiieler 
(1989),  the  wake  component  becomes  important  when  the  density  of  the  slow  stream  is  much 
larger  than  that  of  the  fast  stream.  In  such  instances,  linear  stability  analysis  shows  that  the  early 
development  of  the  layer  is  dominated  by  the  wake  mode  whose  amplification  rate  is  higher  than 
that  of  the  shear  layer  mode.  In  the  renudning  cases,  the  shear  layer  mode  is  dominant,  and  leads 
to  the  familiar  roll-up  of  die  Kelvin-Helmholtz  waves. 

While  the  wake  component  might  be  neglected  in  2D  models,  results  of  Ashurst  &  Meiburg 
(1988)  have  shown  that  the  detail  of  the  vorticity  distribution  plays  an  important  role  in  the 
development  of  3D  instability  modes.  Using  two  vorticity  layers  of  opposite  sign,  they  predicted 
an  asymmetric  spacing  of  the  streamwise  vortex  rods  similar  to  that  experimentally  observed 
(Lasheras  &  Choi,  1988).  However,  this  approach  is  conqilicated  by  the  difficulty  in  specifying 
the  initial  strength  and  separation  of  the  vorticity  layers.  In  fact,  the  initial  separation  of  the 
individual  vorticity  layers  is  not  uniquely  determined  in  that  model,  and  large  separation  distances 
may  lead  to  die  independent  roll-up  of  each  layer,  a  behavior  that  is  not  obtained  in  shear  layers. 

In  our  study,  consideration  of  shear  layers  with  asymmeoic  vorticity  profiles  is  motivated 
by  the  results  of  the  variable-density  layer.  In  particular,  the  numerical  experiment  is  designed  to 
mimic  the  early  effects  of  density  variation  which  were  ^own  to  promote  die  spanwise  voiticity  of 
the  low-density  stream  at  the  expense  of  that  on  the  high-density  side  through  asymmetric  strain. 
This  imitation  is  obtained  by  perturbing  the  second-order  Gaussian  vorticity  profile  to  yield  the 
asymmetric  vorticity  distribution  shown  by  a  broken  line  in  figure  1.  A  single  vorticity  layer  in 
which  the  velocity  increases  mraotonically  &om  one  stream  to  the  other  is  considered.  The  scalar 
is  assumed  passive,  and  its  initial  profile  is  the  saitK  as  diat  used  in  Section  3.1.  We  keep  the  same 
dimensions  and  boundary  conditions,  and  apply  the  same  perturbation  at  the  start  of  the 
confutations.  Since  the  deviation  from  the  symmetric  Gaussian  profile  is  small,  we  expect  similar 
growth  characteristics  of  the  2D  component  of  the  perturbation,  although  its  wavelength  may  not 
correspond  to  diat  of  the  most  amplified  Kelvin-Helmholtz  mode. 

3.3.1  DEFORMATION  OF  MATERIAL  SURFACES 

Perspective  views  of  the  material  surface  initially  lying  in  the  plane  z-0,  shown  in  figure 
21,  exhibit  qualitative  similarity  to  that  of  the  first  case  and  hence  a  detailed  analysis  of  the 
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Lagrangian  motion  is  omitted.  Our  discussion,  will  be  restricted  to  the  features  by  which  the 
asymmetric  layn  may  be  identified.  The  deformation  of  the  material  layers  is  asymmetric  with 
respect  to  the  surface  initially  at  z  =  0,  as  indicated  in  Table  3.  This  development  resembles  that 
observed  in  the  variable-density  case  where  baroclinic  vorticity  generation  produces  a  similar 
asymmetric  strain  field.  Despite  this  similarity,  the  the  two  cases  can  be  distinguished  by  the  fact 
that,  in  the  asymmetric  uniform-density  layer,  the  core  does  not  exhibit  any  convective  motion. 
The  cross-sections  of  the  material  surfaces,  plotted  in  figure  22,  show  that  the  asymmetry  is 
primarily  exhibited  in  the  core  where  the  top  and  bottom  mushrooms  are  distinguishable. 
However,  unlike  the  previous  case,  diere  is  no  indicadtm  that  th^  have  different  strengths. 

Tables 

Normalized  surface  area.  A,  and  number  of  elements,  N,  for  die  individual  material  layers. 


t  =  4.0 

t  =  8.0 

t=12.0 

t  =  16.0 

Material  Layer  1 
Location:  z  =  -132 

A  =  1.03 
N^l.00 

A  =  1.14 
N=120 

A  =  1.72 
N=1.69 

A  =  354 
N  =  4.16 

Material  Layer  2 
Location:  z  s  -0.66 

A  =  1.04 

N  =  1.00 

A  =  124 

N  =  133 

A  =  2.40 

N  =  234 

A  =531 
N  =  855 

Material  Layer  3 
Location:  z  s  0.0 

A  =  1.04 

N  =  1.00 

A  =  1.44 

N  =  1.44 

A  =  2.78 

N  =  3J1 

A  =  5.87 
N  =  9.95 

Material  Layer  4 
Location:  z  s  0.66 

A  =  1.04 
N=1.00 

A  =  122 
N=131 

A  =236 

N  =  2.48 

A  =  5.01 
N  =  735 

Material  Layo:  5 
Location:  z  =  132 

A  =  1.03 
N=1.00 

A  =  1.11 
N=1.15 

A  =  156 

N  =  1.73 

A  =  338 
N  =  3.91 

3.3.2  VORTiaTY  AND  SCALAR  DISTRTOUTION 

Figures  23-26  show  scalar  and  vorticity  contours  plotted  in  the  non-linear  stages  of 
evolution  of  the  flow,  t  =  12J0  and  16,0.  The  spanwise  sections  show  some  but  not  all  the  effects 
associated  with  density  variation  exist  in  this  case.  While  the  convective  motion  could  not  be 
ct^itured,  the  voiticity  and  scalar  distribudons  exhilnt  weak  asymmetry  at  all  qianwise  cuts,  and  the 
rotation  of  the  core  resembles  die  eccentric  rotation  of  an  oval-shiq;>ed  body.  Asymmetric 
entrainment  patterns,  whereby  more  fluid  from  the  bottom  stream  reaches  the  core,  are  also 
established.  At  r  =  76.0,  the  composition  of  the  eddy  slightly  favors  bottom  layer  fluid,  as  the  net 
entrainment  ratios  admit  small  deviations  from  unity,  (SsoTliSso)*  =  1,026,  and  (S'3d)'I(S’3d)* 
-  1.037.  This  corroborates  our  discussion  of  the  motion  of  the  material  surfaces,  whose 
deformation  is  found  more  severe  in  die  bottom  stream.  On  the  other  hand,  the  similarity  between 
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die  asymmeiric  and  varialde^nsity  layers  is  restricted  to  the  fact  that  these  layers  can  be  idoitified 
by  larger  qianwise  vorticity  values. 

The  ginning  of  the  asymmetric  spanwise  eddy  causes  the  generation  of  streamwise  vortex 
rods  of  different  strengths  as  they  approach  the  core.  Unlike  the  variable-densiQ^  case,  the  vortex 
rods  on  the  bottom  side  are  intensified  at  a  lower  rate  in  the  early  stages,  t  <  12.0,  and  this  effect  is 
reversed  at  the  later  stages.  At  r  s  12.0,  the  rado  of  the  circuladon  of  the  two  streamwise  vortices 
is  Ffopirifoi  *  1J13,  and  decreases  to  rtoplFbot  ~  0.982  at  r  =  16.0.  This  asymmetry  occurs 
mainly  as  a  result  of  the  difference  in  shape  and  size  of  the  streamwise  rods  as  the  maximum 
vorticity  values  in  both  remain  close.  A  simplified  model  clarifying  the  origin  of  the  difierence  in 
strengdis  between  the  streamwise  rods  can  be  constructed  by  noting  diat  the  stretching  of  the 
streamwise  vorticity  occurs  along  the  streamwise  boundaries  of  the  domain  where  the  braids  are 
anchored  and  pulled  towards  the  core  of  the  eddy.  In  die  braids,  the  strain  is  weakly  dependent  on 
the  detail  of  the  distribution  within  the  core.  In  fact,  the  strain  in  the  braids  may  be  approximated 
by  concentrating  the  spanwise  vorticity  of  the  core  along  its  center.  Thus,  the  asymmetry  of  the 
vorticity  distribution  of  the  spanwise  core  and  of  its  induced  flow  and  strain  fields,  are  not 
expected  to  play  a  major  role  in  the  production  of  streamwise  vorticity.  The  observed  asymmetry 
becomes  essentially  a  consequence  of  the  eccentric  qnnmng  of  the  spanwise  cores. 

In  the  variable-density  layer,  die  above  description  is  modified  by  the  baroclinic  production 
of  vorticity.  The  streamwise  component  of  die  baroclinic  torque,  r^,  is  small  in  the  braids,  so  that 
the  production  of  streamwise  vorticity  there  is  dominated  by  strain  and  the  contribution  of  the 
baroclinic  torque  to  the  strength  of  the  streamwise  vortices  only  appears  as  a  modulation  of  the 
local  vorticity  values.  On  the  other  hand,  increases  significantly  as  we  approach  the  core. 
Focusing  on  the  later  stages  of  the  flow,  this  observation  is  justified  by  realizing  that  fluid 
acceleration,  which  controls  the  production  of  baroclinic  vorticity,  becomes  dominated  by  the 
convective  accelmtion.  In  turn,  the  convective  acceleration  must  be  larger  in  regions  of  higher 
curvature.  Therefore,  in  the  variable-density  layo:,  as  the  streamwise  vortices  approach  the  core , 
their  strengths  is  affected  by  the  baroclinic  source  tnm  which  is  shown  to  strengthen  (weaken)  the 
streamwise  vorticity  of  the  vortex  rods  lying  on  the  low-density  (high-density)  side. 

Our  numerical  results  have  shown  that  three  mechanisms  can  lead  to  symmetry  loss  in  3D 
tenqwral  shear  layers.  The  first  is  associated  with  the  growth  of  the  translative  instability  which 
leads  to  a  preferential  entrainment  pattern  revving  itself  every  half  spanwise  wavelength  of  the 
instability.  The  second  is  caused  by  an  initially  asymmetric  vorticity  distribution  which  induces 
asymmetric  entrainment  currents  favoring  die  free  stream  having  higher  spanwise  vorticity  values. 
Finally,  baroclinic  vorticity  generation  plays  an  important  role  in  the  formation  of  asymmetric 
structures  by  imparting  a  convective  motion  on  the  cores  and  redistributing  the  vtnticity  within  the 
cores.  These  mechanisms  should  be  distinguished  fiom  similar  effects  occurring  in  2D  and  3D 
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spatially-developing  layers,  in  which  individual  mechanisms  may  be  hard  to  isolate  since 
asymmetric  entrainment  favoring  the  high-speed  stream  is  observed  before  the  transition  to  3D 
motion,  and  the  location  of  vortex  amalgamation,  which  lead  to  the  formation  of  composite 
structures  of  even  higher  complexity,  is  not  easily  predicted  (Dimotakis,  1986). 

The  use  of  asymmetric  vorticity  distribution,  with  a  strong  negative  component,  has  been 
suggested  by  Ashurst  &  Meiburg  (1988)  as  a  means  of  simulating  the  effect  of  the  velocity  ratio. 
Their  computations  were  able  to  C2y)ture  the  asymmetric  spacing  of  the  streamwise  vortices, 
observed  experimentally  by  Lasheras  &  Qioi  (1988).  This  asymmetric  spacing  was  not  observed 
in  our  computations  which  did  not  include  a  negative  initial  vorticity  component.  Thus,  we 
conclude  that  the  presence  of  negative  vorticity  alone  cannot  be  re^nsible  for  this  effect  and  the 
asymmetric  spacing  of  the  streamwise  vortices  is  not  a  property  of  the  "primary"  3D  structure. 
However,  the  development  described  by  Ashurst  &  Meiburg  (1988)  could  occur  during  the 
merging  of  the  distorted  q>anwise  eddies.  In  view  of  the  results  of  Metcalfe  et  al.  (1987),  who 
showed  that  in  the  absence  of  pairing,  both  2D  and  3D  instabilities  tend  to  saturate,  the  suggestion 
that  the  asymmetric  spacing  is  a  product  of  a  higher-order  instability  could  be  justified.  This 
interpretation  is  consistent  witii  the  experimental  findings  of  Bernal  &  Roshko  (1986),  and  with  the 
numerical  simulations  of  Grinstein  et  al.  (1989),  who  observed  that  an  asymmetric  reorganization 
of  the  streamwise  vortices  occurs  after  the  merging  of  the  distorted  spanwise  eddies. 
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5.  CONCLUSIONS 


The  tnmsport  element  method  was  applied  to  study  the  evolution  of  temporal,  doubly- 
periodic,  uniform-  and  variable-density  shear  layers.  The  numerical  schemes  are  Lagrangian  and 
adaptive.  They  are  based  on  tracking  the  vortkity.  scalar  and  scalar  gradients,  discretized  using  a 
finite  number  of  computational  elements.  Three  cases  are  considered:  (i)  a  uniform-density 
vorticity  layer  with  symmetric  vorticity  profile,  (ii)  a  variable-density  vorticity  layer,  and  (iii)  a 
uniform-density  asymmetric  layer,  hi  the  uniform-density  case,  we  take  advantage  of  Kelvin's 
circulation  theorem  in  order  to  save  computational  effort,  while  the  vorticity  transport  equation  is 
used  when  the  beroclinic  source  tom  is  present  Thus,  the  vtnlBX  stretching  is  implicitly  taken  into 
account  in  the  first  case,  while  direct  evaluation  of  tiie  vorticity  stretching  and  baroclinic  source 
terms  is  performed  in  the  second. 

Starting  from  equal  "small-anplitude"  two-  and  three-dimensional  perturbations,  the 
evolution  of  the  vorticity  layer  first  exhibits  a  2D  regime  which  is  characterized  by  the  growtii  and 
roll-up  of  the  Kelvin-Helmholtz  mode.  Following  this  stage,  3D  perturbations  are  rapidly 
amplified.  In  all  cases,  and  consistent  with  previous  results,  two  types  of  3D  instability  are 
observed:  an  instability  in  the  todds  which  leads  to  the  formation  of  streamwise  vortex  rods  and 
scalar  mushroom  structures,  and  a  core  instability  which  causes  an  uneven  deformation  of  the 
spanwise  eddies.  The  instability  in  die  braids  is  a^ociated  with  the  sevoe  stretching  of  the  vortex 
lines  whose  extension  in  the  streamwise  direction  exceeds  the  separation  distance  between 
neighboring  spanwise  eddies.  The  streamwise  vortex  rods  are  continuously  wrapped  around  the 
deformed  spanwise  cores  leading  to  die  intensification  of  the  streamwise  and  spanwise  components 
of  vorticity  and  to  the  generation  of  complex  vortex  structures.  While  the  amplitude  of  instability 
of  the  spanwise  eddies  does  not  reach  such  large  values,  it  is  still  found  to  play  an  important  role  in 
the  evolution  of  die  flowfiekL 

A  detailed  visualization  of  the  vorticity  and  scalar  fields  and  of  the  motion  of  material 
surfaces  was  performed.  The  study  focuses  on  the  manifestation  of  3D  instabilities  in  vorticity 
layers.  The  3D  instalnlity  in  a  vorticity  layer  with  initially  symmetric  vorticity  and  scalar  profiles 
exhibits  asymmetric  vorticity  and  scalar  distributions  at  different  spanwise  locations.  These  forms 
of  asymmetry,  which  reverse  themselves  every  one-half  spanwise  periodicity  length,  are  linked  to 
the  development  of  the  braids  and  translative  instabilities.  Similar  entrainment  currents  are 
observed  in  a  uniform-density  asymmetric  vorticity  layer,  where  the  asymmetry  of  the  flowfield 
leads  to  a  preferential  entrainment  of  irrotational  fluid  from  the  stream  initially  having  higher 
spanwise  vorticity  values.  In  this  case,  the  effea  the  asymmetric  vorticity  distribution  combines 
with  that  of  the  3D  instability,  whose  development  is  not  significanUy  altered  from  the  previous 
case,  resulting  in  a  net  departure  from  a  unity  entraiiunent  ratio.  The  variable-density  layer  is 
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identified  by  a  finite  convective  speed  of  the  eddies  in  the  direction  of  the  high-density  fluid  stream 
and  by  an  asymmetric  entrainment  pattern  favoring  the  low-density  stream.  Baroclinic  torques 
affect  the  development  of  the  3D  nKxles  via  uneven  intensification  or  weakening  of  the  streamwise 
vorticity. 

Special  attention  was  paid  to  the  evolution  of  the  unstable  modes,  and  on  their  roles  in 
reorganizing  the  flowfield.  We  were  thus  content  to  observe  unstable  modes  given  by  a  single 
spanwise  wavelength,  which  was  chosen  close  to  the  most  amplified  mode  of  the  linear  stability 
theory.  We  leave  to  a  subsequent  study,  the  task  of  trying  to  determine  the  characteristics  of  the 
processes  which  lead  to  die  wavelength  selection  of  the  3D  modes.  Moreover,  we  have  restricted 
our  study  to  the  formation  and  maturation  of  the  primary  3D  structures.  We  have  dius  omitted  the 
pairing  interactions  between  these  structures  ami  their  role  in  the  mixing  transition  and  growth  of 
the  layer.  These  confutations  are  currendy  considered. 
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nOURE  CAPTIONS 


Figim  1.  (a)  Schematic  r^nesentation  of  the  shear  layer,  the  coordinate  axes,  and  the  initial 
vordcity  and  scalar  distributuMis,  showing  the  shape  of  the  pertmbatitm,  and  the  initial  location  of 
the  vortex  tubes,  (b)  Vordcity  and  velocity  profiles  for  die  symmetric  and  asymmetric  layers 
discussed  in  Secdons  3.1  and  3.3  respecdvely. 

Figure!.  Three-dimensional  persprodve  view  of  die  s  %  0,  initially  lying  in  the  plane  z  =  0.  The 
plots  were  generated  from  the  point  of  view  of  an  obs^er  loca^  at  (48,24,48).  x-  is  the 
streamwise  direcdon,  y-  the  qianwise  direction,  and  z-  is  the  cross-stream  direction. 

Figure  3.  Three-dimensional  peispecdve  view  of  die  iso-scalar  surface  (a)  initially  lying  in  die 
pli^  z  =  -132,  (b)  embedded  in  die  middle  layer  and  coiiiciding  with  tte  axis  of  the  spanwise 
ctxt,  and  (c)  embedded  in  the  middle  layer  and  located  within  the  bnids.  The  plots  are  generated 
as  in  Figure  1. 

Figure  4.  Intersecdon  of  die  Lagragian  mesh  at  r  s  18.0  with  the  idanes  (a)  y  =  33,  (b)y-  1.6, 
(c)  X  =  6.6,  and  (d)  x  =  2.0.  The  intersecdon  points  are  illustrated  in  terms  of  small  circles  whose 
radius  is  1/8  of  Ae  core  radius  of  the  tran^xirtelanents. 

Figure  5.  Schemadc  illustradon  of  the  evoludon  of  die  scalar  disttibudon  in  the  streamwise  plane 
dividing  the  spanwise  eddy  core. 

Figure  6.  (a)  Contours  of  constant  qianwise  votticity,  a>y,  and  (b)  iso-scalar  contours  plotted  in  tile 
plane  y  =  33. 

Figure  7.  (a)  Contours  of  constant  qianwise  vordcity,  oi^  and  (b)  iso-scalar  contours  jdotied  in  the 
x-z  plane  y  s  5.0. 


Figure  8.  (a)  Contours  of  constant  streamwise  vOTtici^,  Ojc,  and  (b)  iso-scalar  contours  plotted  in 
the  plane  X  s  d.6. 

Figure  9.  (a)  Contours  of  constant  streamwise  vordcity,  otc,  and  (b)  iso-scalar  contour^otted  in 
the  plane  X  =  2.0 . 


Hgure  10.  Evoludon  of  the  total  streamwise  vordci^,  Iloa^  dA,  computed  in  the  streamwise  plane 
located  at  X  =  d.d. 

Hgurell.  Schematic  illustration  of  die  eddy. 

Figure  12.  Evolution  of  the  eddy  size  in  two  dimoisions,  Szd  (v),  and  in  three  dimoisions,  Ssd 

(A),  and  S'iD(*). 

Figure  13.  Three-dimensional  perspective  view  of  die  surface  p' s  0,  initially  lying  in  the  plane  z 
=  0  for  the  case  of  a  variable-densi^  shear  layer.  The  plots  were  generated  firom  the  point  of  view 
of  an  observer  located  at  (48,24,48). 

Figure  14.  Intersection  of  the  Lagrangian  mesh  at  r  s  16.0  widi  the  planes  (a)  y  =  33,  (b)  y  =  1.6, 
(c)  X  *  9J,  and  (d)x  =  2.0. 


Figure  15.  Contours  of  spanwisc  vorticity,  tOy,  (left)  and  density  (right)  plotted  in  the  plane  y  = 
33. 

Figure  16.  Contours  of  spanwise  vorticity,  (i>y,  (left)  and  density  (right)  plotted  in  the  plane  y  = 
5.0  . 

Figure  17.  Contours  of  streamwise  vorticity,  tUr,  (left)  and  density  (right)  plotted  at  r  =  12.0,  and 
16.0,  in  the  planes  x  =  8.6,  and  jc  =  9  J  respectively. 

Figure  18.  Contours  of  streamwise  vorticity,  0)x,  (left)  and  density  (right)  plotted  in  the  plane  x  = 

2.0. 

Figure  19.  Contours  of  the  spanwise  component  of  the  baroclinic  torque,  Ty,  plotted  in  the  plane 

(a)  y  =  i  J  and  (b)  y  =  5.0. 

Figure  20.  Contours  of  the  streamwise  component  of  the  baroclinic  torque,  Tx,  plotted  at  r  =  8.0 
and  16.0,  in  y-z  planes  given  by  (a)  x  =  8.0  and  93;  and  (b)  jc  =  2.0. 

Figure  21.  Thiee*dimensional  perspective  view  of  the  surface  s  ^0,  initially  lying  in  the  plane  z  = 
0  fOT  the  shear  layer  with  an  initially  asymmetric  vorticity.  The  plots  were  general^  from  the  point 
of  view  of  an  observer  located  at  (48,24,48). 

Figure  22.  Intersection  of  the  Lagrangian  mesh  at  r  s  16.0  witii  the  planes  defined  by  (a)  y  s  33, 

(b)  y  =  1.6,  (c)  X  =  6.6,  and  (d)  x  =  2.0. 

Figure  23.  Spanwise  vorticity,  <Uy,  anti  scalar  contours  (right)  plotted  in  tiie  plane  y  s  5  J . 
Figure  24.  Spanwise  vorticity,  (Oy,  (left)  and  scalar  contours  (right)  plotted  in  the  plane  y  =:  5.O.. 
Figure  25.  Streamwise  vorticity,  (%,  Oeft)  and  scalar  contours  (right)  plotted  in  the  plane  x  =  6.6. 
Figure  26.  Streamwise  vorticity,  cojc,  (left)  and  scalar  contours  (right)  plotted  in  the  plane  x  =  2.0. 
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ABSTRACT 


The  three-dimensional  transport  element  method  is  extended  to  solve  the 
conservation  equations  for  reacting  flow.  The  numerical  scheme  belongs  to  an 
adaptive,  Lagrangian  class  of  field  methods  in  which  computational  effort  is 
concentrated  in  zones  of  finite  vorticity  and  chemical  reaction.  The  method 
is  based  on  accurate  discretization  of  vorticity  and  species  concentrations 
among  a  number  of  spherically-symmetric  transport  elements,  and  the  advection 
of  these  elements  along  particle  trajectories.  We  use  the  low  Mach  number 
approximation  of  combustion  in  open  domains,  and  restrict  our  attention  to  the 
case  of  diffusion  flames  with  no  heat  release.  A  single-step,  second-order, 
infinite-rate  kinetics  chemical  reaction  model  is  employed. 

The  scheme  is  applied  to  study  the  effect  of  flow-induced  instabilities 
on  the  reaction  zone  and  product  distribution  in  a  temporal  shear  layer. 
Results  are  obtained  in  the  high  Peclet  number  regime  for  a  wide  range  of 
Damkohler  numbers.  Changes  in  the  reaction  field  are  related  to  either  the 
entrainment  or  the  strain  field  associated  with  the  saturation  of  the 
instabilities.  With  increasing  Damkohler  number,  the  structure  of  the 
reaction  region  changes  from  a  distributed  zone  embedded  within  spanwise  and 
streamwise  vortices,  to  a  thin  sheet  surrounding  their  cores.  However,  the 
product  concentration  always  exhibits  strong  similarity  to  the  vorticity 
distribution,  realizing  its  highest  values  in  zones  of  high  vorticity  and 
falling  rapidly  in  regions  where  the  vorticity  magnitude  is  small.  Variation 
of  the  Peclet  number  yields  minor  changes  in  the  product  distribution  and  in 
the  reaction  zone  structure,  but  strongly  affects  product  formation  rates. 


I  INTRODUCTION 


The  subject  of  this  work  is  the  construction  of  Lagrangian  field  methods 
for  the  numerical  simulation  of  three-dimensional,  high  Reynolds  nvimber, 
reacting  shear  flow.  In  this  fluid  flow  regime,  vorticity  is  confined  to  a 
small  well-defined  fraction  of  the  flow  field  making  vortex  methods,  in  which 
computational  elements  are  used  to  cover  the  support  of  vorticity,  a  natural 
candidate  for  the  solution  of  the  momentxjm  equation.  In  previous  effort, a 
modified  vortex  element  scheme  was  constructed,  atnd  applied  to  study  the 
motion  of  \inconfined,  inviscid  vortex  rings.  The  study  revealed  the  need  for 

(1)  a  careful  representation  of  the  vorticity  field  at  the  initial  state,  and 

(2)  an  accurate  inplementation  of  numerical  schemes  as  time  evolves.  These 
requirements  are  especially  critical  in  three-dimensional  flows  where  several 
forms  of  rapidly-growing  instabilities  of  complex  shapes  are  present.  To 
satisfy  these  requirements,  and  to  accommodate  the  solution  of  scalar 
conservation  equations,  an  adaptive  class  of  Lagrangian  field  methods,  called 
transport  element  methods,  was  constructed.  Numerical  study  of  several  forms 
of  these  methods,  which  focused  on  the  solution  of  the  inviscid  scalar 
transport  equation,  was  successfully  conducted  in  Refs.  4,5. 

While  the  ultimate  goal  of  this  work  consists  of  the  construction  of 
numerical  methods  for  the  solution  of  the  compressible,  turbulent  combustion 
equations,  in  this  stiidy,  we  limit  our  attention  to  the  low  Mach  nutdjer, 
isothermal  flow  limit.  In  the  presence  of  an  exothermic  reaction,  analysis  of 
burning  rates  is  complicated  by  the  expansion  field  and  the  baroclinic 
vorticity  generation  spawned  by  the  heat  release,  aixJ  the  burning  rates 
exhibit  higher  sensitivity  to  the  strain  field.  In  a  three-dimensional  flow, 
these  mechanisms  are  further  compounded  by  vorticity  stretching  and  individual 
effects  become  hard  to  isolate.  By  virtue  of  our  simplifying  assumptions. 


computed  results  will  allow  us  to  focus  on  studying  the  deformation  of  the 
diffusion  flame  by  a  known  flow  field. 

The  numerical  scheme  is  obtained  by  modifying  the  transport  element 
method  discussed  in  Ref.  5.  It  is  based  on  the  discretization  of  both 
vorticity  and  species  concentrations  into  a  number  of  transport  elements  of 
finite,  spherically-symmetric  overlapping  cores.  Transport  elements  are 
distributed  over  entire  material  surfaces  whose  motion  is  tracked  in  a 
Lagrangian  frame  of  reference.  Vorticity  associated  with  the  elements  changes 
by  stretching  and  tilting,  while  species  concentrations  are  updated  by 
numerically  accounting  for  diffusion-reaction  terms.  The  rotational  velocity 
field  is  confuted  by  a  discrete,  desingular i zed  Biot-Savart  convolution  over 
the  field  of  the  elements.  Boundary  conditions  are  satisfied  by  adding  the 
contribution  of  the  appropriate  image  system  of  the  transport  elements. 

The  method  is  applied  to  study  the  evolution  of  chemically- reacting 
temporal  shear  layers.  The  flow  field  is  selected  because  of  (1)  its 
practical  importance  in  a  large  class  of  combustion  problems,  and  (2)  the 
existence  of  experimental®”^®  and  analytical^^”^^  results  which  focus  on  the 
three-dimensional  behavior  of  the  flow.  However,  numerical  studies  of  scalar 
entrainment,  mixing  and  chemical  reaction  in  three-dimensional  layers  remain 
scarce.  ”  Thus,  the  method  is  used  to  study  the  effect  of  vorticity- 
induced  entrainment  and  strain  on  the  evolution  of  diffusion  flames.  Computed 
results  are  obtained  for  a  wide  range  of  Oamkohler  numbers,  in  order  to  study 
its  effect  on  the  structure  of  the  reaction  zone  and  product  distribution. 

II  FORMULATION  AND  NUMERICAL  SCHEME 


II.l  FORMULATION  AND  GOVERNING  EC3UATI0NS 


We  consider  a  three-dimensional,  chemically-reacting  flow  in  the  limits 
of  infinite  Reynolds  number,  vanishingly  small  Mach  number  eind  heat  release. 
The  reaction  is  described  by  a  single-step,  second-order,  infinite-rate 
kinetics,  F  (fuel)  +  0  (oxidizer)  -♦  P  (products).  The  reactants  and  products 
are  assumed  to  behave  as  perfect  gases  with  equal  molecular  weights,  equal  eind 
constant  physical  properties.  Under  these  assumptions,  the  vorticity  form  of 
the  momentum  equation,  the  continuity,  and  species  conservation  equations, 
normalized  with  respect  to  the  appropriate  combination  of  velocity,  U^, 
length,  L^,  eind  density,  p^,  scales,  reduce  to: 


Du 

Dt 


7*u  -  0 


(1) 

(2) 


Ds^ 

Dt 


1 

Pe  Le 


(3) 


In  Eqs.  (l)-(3),  u  ■  (u,v,w)  denotes  the  velocity  vector  in  a  right-handed 
Cartesian  system  x  -  (x,y,2),  «•  Txu  is  vorticity,  t  is  time,  D/Dt  -  3/3t  + 
u*7  is  the  material  derivative,  7  ■  0/3x,3/3y, 3/3z)  is  the  gradient  operator, 
7^  is  the  Laplace  operator,  Pe  -  U^iaL^)  is  the  Peclet  number,  Le  •  a/fi  is 
the  Lewis  nvanber,  vrtiile  a  and  0  denote,  respectively,  the  thermal  and  mass 

diffusivities.  The  model  requires  the  transport  of  two  scalars,  s^  i-1,2, 

1  2 
where  s  ■  Cq  denotes  the  oxidizer  concentration,  while  s  ■  Cp  represents  the 

fuel  concentration.  The  product  distribution  is  given  by  Cp  ■  1  -  Cq  -  Cp, 

and  chemical  production  term  W  «  Da  c^  Cp  where  Da  is  the  Damkohler  number. 

Equations  (l)-(3)  are  obtained  by  simplifying  the  low  Mach  nundaer  equations 

for  condoustion  for  vanishingly  small  heat  release.  The  assumption  of  no  heat 

release  implies  that  both  tenperature  and  density  remain  constant,  so  that 

both  baroclinic  vorticity  generation  and  volumetric  expansion  effects  cancel 

in  the  vorticity  transport  equation.  The  model  ignores  viscous  effects,  since 


viscosity  only  affects  the  motion  of  the  vorticity  field  well  beyond  tr.e 
mixing  transition,  a  mechanism  '//hich  is  not  in  the  scope  of  the  present 
computations.  Mass  diffusion  effects  are  retained,  because  they  govern  the 
mixing  and  chemical  reaction  processes. 


II. 2  NUMERICAL  SCHEME 

The  numerical  scheme  used  in  the  solution  of  the  governing  equations  is 
obtained  by  extension  of  the  methods  analyzed  in  previous  work,^  and  hence  is 
summarized  in  the  following.  Its  construction  starts  with  the  discretization 
of  the  vorticity  field,  (a(x) ,  on  a  three-dimensional  mesh: 

N 

(D(x,0)  -  £  w.  (0)  dV.  fj.(x-X.  )  (4) 

i^l  1  15  1 

where  N  is  the  total  nvunber  of  vortex  elements,  Xj;,  dv^,  and  denote  the 
center,  volume,  and  vorticity  of  element  i,  respectively.  The  vorticity 
associated  with  each  element  is  smoothed  in  a  small  neighborhood  of  X^ 
according  to  a  spherical  core  function,  f^,  with  core  radius,  i,  where  fj(r)  » 
(1/8^)  f(r/8).  A  third-order  Gaussian  core  function, 

f(c)  -  e"*^  (5) 

18^19 

is  adopted.  Note  that  f  decays  rapidly  for  r  >  8  so  that  8  represents 

the  radius  of  the  sphere  where  vorticity  is  concentrated. 

The  vortex  elements  are  initialized  by  employing  entire  material  surfaces 
distributed  within  the  support  of  vorticity.  The  surfaces  are  discretized 
into  rectangular  "transport"  elements  which  generalize  the  notion  of  the 
vortex  vector  elements.  In  this  representation,  a  transport  element  is 
specified  by  the  vector  where  the  Lagrangian  coordinates 

X^,  and  are  used  to  define  a  "rectangular”  area  around  its  center.  The 
transport  elements  are  selected  such  that  the  sides  of  adjacent  rectangles 


coincide  to  form  a  mesh  which  continuously  describes  an  entire  material 
surface.  This  construction  resembles,  but  differs  from  that  proposed  by 
Agishtein  &  Migdal  who  used  triangular  elements  to  describe  a  singular  vortex 


sheet.' 


Within  each  element,  we  use  a  finite  element  description  of  the  surface 

..  1.*....^..  ^  ^  ^  .  21  ml. ..  — £  ..I . _ 


using  linear  interpolation  functions.' 


The  motion  of  the  centers  of  the 


elements,  X^,  is  approximated  by: 

Xi(t)  -  \  (xj(t))  +  ;^(t)  +  x^(t)  +  xj(t)) 


where  the  ;^(t),  j  -  1,2, 3, 4,  denote  the  instantaneous  coordinate  of  the 

material  particles  whose  Lagrangian  coordinates  coincide  with  the  vertices  of 

the  initial  mesh.  This  representation  accommodates  the  integration  of  the 

vorticity  transport  equation,  v^ich  is  necessary  in  general  flow  situations,^ 

but  retains  the  advantages  of  vortex  element  schemes  which  relate  the 

22  26 

evolution  of  the  vorticity  field  to  the  deformation  of  vortex  lines. 

Under  the  assumptions  made  above,  the  latter  approach  is  preferred,  since  it 
results  in  considerable  computational  savings  as  direct  evaluation  of  the 
gradient  of  the  flow  map  is  avoided.  This  is  done  by  associating  with  each 
element  a  circulation  and  requiring  that  the  pair  of  opposing  sides  of  the 
rectangles,  (X^,X^)  and  align  with  the  local  vorticity  vector.  The 

quantity  «^(t)dV^  is  thus  replaced  by  r^5Xj^(t),  where  the  evolution  of  the 
material  length  8Xj^(t)  is  given  by: 

&Xi(t)  -^(X^(t)  +  )^(t)  -  xj{t)  -  X?(t))  (7) 


According  to  Kelvin's  theorem,  F^^ 
while  Helmholtz'  theorem  is  used  to 


8Xj^(t)  as  follows: 

»4  (t)  »  - 

^  l«il 


8Xi(t) 


remains  constant  along  a  particle  path, 
relate  the  evolution  of  <a^{t)  to  that  of 


(8) 
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The  velocity  field,  induced  by  the  vorticity  distribution  in  an  infinite 
domain,  is  given  by  the  discrete,  desingularized  Biot-Savart  law: 


where  K(r)  «  4n  f(r')  r'^  dr'  -  l-exp(-r^),  and  r^  »  |x-X^|.  The  velocity 
field  thus  obtained  is  used  in  conjunction  with  a  second-order  predictor- 
corrector  scheme  to  update  the  particle  positions  X^(t). 

Analysis  of  particle  methods,  2uid  nimerical  evidence^  indicate  that 
severe  deformation  of  the  Lagrangian  mesh  under  the  action  of  the  strain  field 
causes  a  deterioration  in  the  discretization  accuracy.  To  overcome  this 
difficulty,  the  remeshing  algorithm  suggested  in  Ref.  5  is  utilized.  It 
effectively  amounts  to  the  redistribution  of  the  elements  along  a  vortex  tube 
and  the  splitting  of  the  vortex  tubes  along  material  surfaces,  whenever  the 
length  of  the  elements  exceeds  the  core  radius.  This  allows  us  to  capture 
severe  distortions  of  the  flow  without  losing  accuracy. 

The  numerical  solution  of  Eq.  (3)  is  performed  in  a  similar  way  as  that 
of  the  vorticity  field.  The  concentration  fields  are  discretized  among  the 
transport  elements,  which  now  carry,  along  with  vorticity,  discrete  local 
values  of  s^  and  s^.  We  let: 


N 

s(x,t)  •  s^i.t)  dv^  fj(x-Xj^(t))  (10) 


Next,  we  adopt  the  discrete  approximation  of  the  Laplacian  operator: 
-  B(5j)  -  ^  (s^  -  Sj)  95(X,(t)  -  x,(t)) 


vdiere 

^  ^  9(-^) 

0 


(12), 
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and  f  is  the  core  smoothing  function  previously  defined.*^  ~  This  scheme  is 
preferred  over  random  walk  methods, whose  implementation  is  cumbersome  in 
conjunction  with  the  transport  of  connected  surfaces,  and  over  core  spreading 
techniques, which  no  longer  yield  simple  evolution  equations  for  the  core 
radii.  The  evolution  of  the  scalar  concentration  fields  is  found  by 
numerically  integrating: 

ds^  1  i 

a?  ■ 

using  the  same  method  employed  in  tracking  the  particle  positions. 

Ill  FLOW  GEOMETRY  AND  INITIAL  CONDITIONS 

We  assvime  a  vorticity  layer  of  finite  thickness,  periodic  in  its 
streamwise  x-direction  and  spanwise  y-direction,  and  unconfined  in  the  cross¬ 
stream  z-direction.  The  thickness  of  the  layer  is  expressed  by  2a  where  a  is 
the  standard  deviation  of  the  second-order  Gaussian  curve  which  describes  the 
physical  vorticity  distribution  within  the  layer  at  t»  0.  Letting  2(x)  denote 
this  initial  condition,  we  have:  2y(x)  -  2/(a  /n)exp(-z  /a  ) ,  2^(x)  ■  “ 

0.  The  corresponding  velocity  distribution  U(x)  at  t  -  0  is  given  by:  V(x)  ■ 
W(x)  ■  0,  and  U(x)  »  erf(z/<j),  erf  being  the  error  function.  The  vorticity 
layer  thus  admits  a  streamwise  velocity  difference  au  ■  2.  The  periodicity 
length  in  the  x-direction,  »  13.2  a,  corresponds  to  the  wavelength  of  the 

two-dimensional  most  unstable  roode,^^  vhile  X^  -  X^2  is  selected  close  to  the 

12 

wavelength  of  the  most  amplified  three-dimensional  mode.  The  initial  fuel 
and  oxidizer  concentrations  follow  error-function  type  profiles,  CQ(t«0)  - 
erf(z/a),  Cp(t-O)  -  1-erf (z/o),  so  that  Cp(t-O)  »  0.  Thus,  the  top  stream  and 
bottom  streams  consist  respectively  of  oxidizer  and  fuel,  and  the  initial 


10 


thickness  of  the  vorticity  layer  matches  the  "mixed"  zone,  or  the  zone  of 
finite  chemical  activity. 

The  layer  is  initially  discretized  among  elements  distributed  on  a  grid 
of  20x14x5  points  along  the  x-,  y-,  eind  z-directions  respectively. 
Furthermore,  a  is  chosen  as  a  reference  length  scale  so  that  Ax»Az*0.66,  and 
Ay-0.471.  The  top  stream  velocity,  U(z-»  *),  is  selected  as  a  reference 
velocity  scale,  the  time  step  At  -  0.1,  and  &  »  0,89.  The  circulations  of  the 
elements  are  found  by  matching  discretized  and  assumed  vorticity  values  at  the 
centers  of  the  elements,  while  the  concentration  fields,  s^  «md  s^^,  are 
initialized  by  associating  with  the  transport  elements,  the  values  realized  by 
the  corresponding  "physical"  distribution  at  their  centers.  Normal,  periodic 
boundary  conditions  are  satisfied  by  adding  the  contribution  of  the  image 
system  of  the  trzmsport  elements,  using  the  procedure  detailed  in  Ref.  5. 

The  layer  is  perturbed  at  t  ■  0  using  two  sinewaves  of  the  same  amplitude 
in  the  streanwi.se  and  spanwise  directions,  as  expressed  by  the  transformation 
Zi  ->  Zi  +  esin(2nx^/X^)  +  csin(2ny^/Xy) ,  e  -  0.02  X.  Computations  of  the  non¬ 
reacting  flow  field  are  carried  until  t  -  18.0  to  observe  the  growth  of  two 
and  three-dimensional  instabilities,  vrfiile  those  of  the  reacting  flow  field 
are  stopped  at  t  -  16,0,  i.e.  when  mature  vortical  structures  are  formed. 
Results  are  obtained  for  a  total  of  21  cases,  assuming  infinite  Reynolds 
number,  unit  Lewis  number  and  varying  the  Damkohler  and  Peclet  numbers.  We 
consider  seven  values  of  the  Damkohler  number.  Da  -  0.1,  0.2,  0.4,  1.,  2.5, 
5.0,  and  10.0,  and  three  values  of  the  Peclet  number,  Pe  -  250,  500,  and  1000. 
In  the  following,  we  start  with  a  svonmary  of  the  evolution  of  the  flow  and 
vorticity  fields,  and  then  discuss  the  influence  of  the  vortical  structures  on 
the  shape  of  the  reaction  zone  and  the  product  concentration  distribution. 


IV  FLOW  FIELD  EVOLUTION 


The  development  of  the  shear  layer,  as  represented  by  the  evolution  of 

the  material  surfaces  where  the  vorticity  in  non-zero  (the  flow  is  inviscid), 

and  the  vorticity  field  on  different  spanwise  and  streamwise  locations  are 

excimined.  Similar  flow  configurations  were  analyzed  by  the  present  authors  in 

F-f.  5,  studied  by  Ashurst  and  Meiburg  who  focused  on  the  dynamics  of  vortex 
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filaments  in  both  symmetric  amd  asymmetric  layers,  and  by  Metcalfe  et  al. 
who  used  a  spectral  scheme  to  examine  the  evolution  of  individual  instability 
modes. Hence,  the  main  features  of  the  flow  are  briefly  discussed  in  the 
following. 

Figure  1  aepicts  three-dimensional  perspective  viev;s  of  the  material 
surface  initially  located  at  z  -  0,  at  t  -  12.0  amd  16.0.  This  represents  the 
middle  surface  within  the  shear  layer  where  most  of  the  vorticity  is 
concentrated.  The  plots  are  generated  from  the  point  of  view  of  an  observer 
located  at  (48,24,48). 

At  t  «  12.0,  the  rollup  of  the  vorticity  layer  produces  a  well-defined 

spanwise  eddy  core.  The  saturation  of  the  instability  is  associated  with  the 

redistribution  of  the  vorticity  field  amd  results  in  the  creation  of  a  core 

where  most  of  the  spanwise  vorticity  is  concentrated,  aind  braids  which  are 

constant’’/  strained  under  the  influence  of  the  spanwise  vortex  cores.  The 

amplitude  of  the  spanwise  perturbation  is  significantly  amplified  along  the 
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core,  an  indication  of  the  evolution  of  the  translative  instability.  This 
non-uniform  axial  displacement  of  the  core  is  accompanied  by  an  out-of  phase 
deformation  of  the  braids  region  under  the  influence  of  the  generated 
streamwise  vorticity. 

At  t  «16.0,  the  motion  becones  highly  three-dimensional  as  differences 
along  various  spanwise  stations  become  in^rtant  and  depart  from  the  shape  of 
the  imposed  sinewave  perturbation.  The  stretching  of  the  braids,  which  are 


"anchored"  along  the  boundaries  of  the  domain  and  pulled  towards  the  core, 
leads  to  the  generation  and  intensification  of  streamwise  vorticity.  This 
results  in  the  saturation  of  the  streamwise  vorticity  into  vortex  rods  which 
extend  throughout  the  braids  and  are  wrapped  around  the  spanwise  core.  The 
presence  of  streaunwise  vortex  rods  is  inferred  from  the  spinning  of  material 
surfaces  about  streaimwise  axes  which  are  located  at  the  spanwise  mid-section 
and  boundaries  of  the  domain. 

The  Lagrangian  description  of  the  flow  is  completed  by  considering  two- 
dimensional  cross-sections  of  the  computational  surfaces.  Figure  2  shows 
cross-sections  through  all  the  material  surfaces  at  t  -  18.0.  We  consider 
spanwise  sections  along  two-dimensional  planes  specified  by  (a)  y  -  3.3;  and 
(b)  y  -  1.6,  and  streamwise  sections  through  the  core  and  braid  regions,  in 
the  planes  located  at  (c)  x  -  6.6;  and  (d)  x  »  2.0.  Circles  are  drawn  to  mark 
the  intersection  points  with  the  transport  elements.  The  radius  of  the 
circles  is  smaller  than  the  core  radii  of  the  smoothing  functions. 

The  spanwise  cross-sections  illustrate  the  effect  of  the  translative 
instability  which  causes  a  non-uniform  axial  deformation  of  the  spanwise 
vorticity  core.  The  eddy  core  eddy  is  pushed  upwards  and  in  the  direction  of 
the  top  stream  in  the  "left"  half  of  the  domain,  0  <  y  <  3.3,  while  it  suffers 
an  antisymmetric  deformation  in  the  other  half.  The  core  region  at  most 
sp6mwise  cross-sections  loses  its  symmetry,  as  computational  elements  migrate 
in  the  direction  opposite  to  that  of  the  core  translation.  The  intersection 
of  material  surfaces  with  the  pleme  y  -  3.3,  shows  that  the  braids  thicken 
significantly  at  this  spanwise  location  and  that  they  entrain  irrotational 
fluid  from  both  free  streams.  This  can  be  verified  by  simultaneously 
examining  Fig.  2c  which  shows  the  wrapping  of  the  mushrooms  around  the  core  of 
the  eddy.  The  entrainment  of  mushrooms,  which  are  generated  in  the  braids, 
towards  the  core  region  leads  to  the  formation  of  a  double  structure.  Dark 
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areas  ’^ich  appear  along  the  core  region  correspond  to  the  intersection  of  the 
streantwise  vortex  rods  with  the  plane  of  the  figure. 

Figures  3-6  show  constant  spanwise  vorticity  contours,  plotted  in  two 
spanwise  sections,  y  -  3.3  and  1.6,  and  streantwise  vorticity  contours  in  the 
streantwise  sections  x  =  6.6  and  2.0,  respectively.  Contours  are  generated  at 
times  t  *  12.0  and  16.0.  At  earlier  times,  the  vorticity  field  exhibits  small 
deviations  from  that  obtained  in  a  two-dimensional  flow,  and  hence  is  not 
depicted.  The  growth  of  three-dimensional  perturbations  is  suppressed  during 
the  initial  stages  of  the  Kelvin-Helroholtz  instadiility  and  the  spanwise 
vorticity  remains  essentially  uniform  across  the  layer Weak 
streainwise  structures  are  generated  by  local  tilting  of  the  vortex  lines  into 
the  streantwise  direction,  but  do  not  significantly  affect  the  evolution  of  the 
flow.  This  mechanism  leads  to  the  creation  of  zones  of  alternating  streantwise 
vorticity  whose  locations  and  signs  follow  the  shape  of  sinewave  perturbation. 

Following  the  rollup  of  the  spanwise  vorticity  the  flow  field  undergoes 
a  rapid  transition  to  three-dimensional  motion.  In  the  braids  region,  a 
violent  increase  in  the  amount  of  streantwise  vorticity  is  observed.  The 
streantwise  structure  almost  doubles  in  strength  in  half  the  time  span  it  took 
to  generate  it.  This  behavior  is  expected  since  the  creation  of  the  large 
spanwise  eddy  results  in  larger  strain  rates  in  the  neighborhood  of  the 
stagnation  "lines"  which  anchor  the  braids. 

In  the  meantime,  the  streantwise  vorticity  distribution  along  the  core 
becomes  distinguished  by  the  presence  of  a  top  emd  a  bottom  row  of  counter¬ 
rotating  streantwise  vortices.  These  rows  are  separated  by  a  third,  middle 
row,  which  is  generated  as  a  result  of  the  growth  of  perturbations  along  the 
core  itself  by  the  mechanism  of  the  translative  instability,  and  is  180°  out 
of  phase  with  the  other  two.^^  The  growth  of  perturbations  along  the  core 
leads  to  a  more  dramatic  increase  of  the  streamwise  vorticity,  as  the  middle 
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row  of  vorticeo  becomes  stronger  than  the  top  and  bottom  rows.  The  generation 
of  strong  strearawise  structures  in  the  core  region  is  not  associated  with 
large  amplitude  deformation  in  the  streamwise  direction.  This  may  be 
explained  by  noting  that  the  rollup  of  the  layer,  which  precedes  the  three- 
dimensional  motion,  forces  the  migration  of  the  spanwise  vorticity  from  the 
thinning  braids  into  the  core  region.  Also  note  that  the  intensification  of 
the  streamwise  vortices,  entrained  from  the  braids  into  the  cores  (the  top  and 
bottom  rows),  lags  that  observed  in  the  braids.  Thus,  streamwise  vorticity 
associated  with,  the  vortex  rods  is  generated  in  the  braids  and  then  strained 
towards  the  cores. 

The  transition  to  three-dimensional  motion,  in  its  initial  stage,  does 
not  lead  to  significauit  qualitative  modification  of  the  structure  of  the  layer 
as  perceived  in  two-dimensional  spauiwise  cuts.  The  spanwise  vorticity 
distribution  still  bears  a  striking  resemblance  to  that  observed  in  a  a  two- 
dimensional  flow.  However,  we  note  two  minor  differences:  (i)  the  vortex 
core,  while  essentially  governed  by  the  dynamics  of  the  primary  instability, 
suffers  a  reduction  of  its  cross-section;  and  (ii)  the  vorticity  distribution 
loses  its  symmetry  in  the  pleme  located  at  y  -  1.6,  but  remains  symmetric  at 
the  spanwise  midsection  of  the  domain,  y  -  3.3.  The  first  effect  is  a 
consequence  of  the  growth  of  the  three-dimensional  perturbation  on  the 
spanwise  core  which  necessitates  a  stretch  component  along  its  axis.  The  loss 
of  symmetry,  which  becomes  more  obvious  in  the  later  stages  of  development  of 
the  three-dimensional  instability,  is  discussed  below. 

The  non-linear  stages  of  evolution  of  the  three-dimensional  instedDility 
are  examined  first  in  the  streamwise  planes  of  the  shear  layer.  Figs.  5  and  6. 
The  stretching  of  the  vorticity  lines  along  the  braids  leads  to  the  maturation 
of  the  elongated  vortices  into  round  concentrated  cores. The  flow  field 
induced  by  the  vortex  rods  causes  further  deformation  of  the  material  surfaces 
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resembles,  cut  is  different  from,  the  two-dimensional  pairing  of  vortices. 
This  type  of  distribution  should  be  contrasted  with  that  obtained  in  unstable 
vortex  rings  which  displays  similar  features  at  certain  azimuthal  cross- 
sections.^ 


V  REACTING  SHEAR  LAYER 

We  turn  our  attention  to  the  study  of  diffusion  flames  in  three- 
dimensional  vorticity  layers.  Since  heat  release  effects  are  neglected, 
temperature  and  density  remain  constant,  and  the  reaction  does  not  affect  the 
flow  field.  Results  are  used  to  analyze  the  effect  of  vorticity-inJuced 
entrainment  and  strain  on  the  reaction  zone,  product  distribution,  emd  burning 
rates.  We  start  with  a  reacting  layer  at  low  Damkohler  number  and  low 
diffusion.  Da  «  0.1  and  Pe  *  1000,  then  at  high  Damkohler  number  and  low 
diffusion.  Da  «  10  and  Pe  ■  1000. 

Results  are  shown  in  terms  of  shaded  contours  of  product  concentration, 
and  product  formation  rate,  W.  Product  concentration  is  shown  in  terms  of  six 
different  shades  of  gray,  mapped  to  equal  size  intervals  between  0  and  1.  The 
product  formation  rates  are  first  normalized  by  the  maximum  value  achieved  in 
the  plane  \diere  the  figure  is  generated,  before  the  contours  are  drawn.  In 
these  plots,  dark  areas  highlight  the  zone  of  highest  chemical  activity,  and 
must  not  be  used  as  indication  of  the  actual  formation  rates. 

V.l  REACTING  LAYER  AT  LOW  DAMKOHLER  NUMBER 

Figures  7-10  show  product  concentration  and  reaction  rate  contours  in  the 
planes  y  >  3.3,  y  •  1.6,  x  *  6.6,  and  x  2.0,  respectively,  generated  at  t  « 
12.0  and  16.0.  In  each  figure,  product  concentration  and  reaction  rate 
contours  are  shown  side-by-side.  At  earlier  stages,  t  <  8.0,  examination  of 


the  reaction  field  (not  shown  here)  exhibits  small  differences  from  those  cne 
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obtains  in  a  two-dimensional  field.  The  products  of  reaction  and  the 
zone  of  highest  chemical  activity  surround  the  region  of  highest  vorticity. 
Three-dimensional  effects  are  weak  in  this  stage,  and  can  only  be  detected 
along  the  streeunwise  sections,  which  show  very  small  spanwise  undulations  of 
the  layer. 

With  the  intensification  of  the  streantwise  vortices  in  the  braids  and  the 
amplification  of  the  tramslative  instability  in  the  cores,  significant 
deviations  from  two-dimensional  behavior  are  observed.  At  y  -  3.3,  the 
streamwise  vortex  rods  which  start  to  form  at  the  boundaries  of  the  domain 
enlarge  the  product  concentration  layer  at  this  location.  As  shown  in  Fig. 
7b,  enlargement  of  the  layer  is  a  consequence  of  burning  enhancement, 
associated  with  the  rollup  of  the  surfaces  around  the  streamwise  vortices. 
This  is  contrary  to  the  result  of  two-dimensional  simulations  in  which  the 
strain  field  of  the  spanwise  vortex  cores  leads  to  continuous  thinning  of  the 
product  layer  embedded  within  (or  surrounding)  the  braids.^® 

However,  the  thinning  of  •  the  braids  and  the  entrainment  of  products 
towards  the  spanwise  vortex  core  are  still  observed  in  the  spanwise  section  y 
-1.6,  which  does  not  intersect  any  streamwise  vortex  rod.  Meanwhile,  the 
products  of  reaction  undergo  a  displacement  similar  to  that  of  the  vorticity 
core  (Fig.  8a).  The  translation  of  the  core  towards  the  top  stream  and  the 
migration  of  its  geometric  center  in  the  opposite  direction  destroy  the 
symmetry  of  the  product  concentration.  The  asynmetric  mixing  patterns 
associated  with  this  mechanism  also  affect  the  reaction  rate  distribution 
(Fig.  8b),  since  the  reaction  rate,  which  depends  strongly  on  the  composition 
of  the  reacting  mixture,  drops  rapidly  if  either  oxidizer  or  fuel  become 
deficient.  This  description  persists  survives  the  restructuring  of  the 
vorticity  core  by  the  maturation  of  three-dimensional  instabilities,  so  that 
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the  reaction  rate  and  product  concentration  distributions  at  this  section  are 
in  qualitative  agreement  with  the  results  of  two-dimensional  simulations."^^ 

The  difference  between  the  two  sections  is  directly  related  to  the 
streamwise  vorticity.  The  presence  of  streamwise  vort=x  rods  leads  to  mixing 
and  combustion  enhemcement,  a  mechanism  which  resists  the  "negative"  effects 
of  the  underlying  two-dimensional  strain  field.  Not  only  do  the  streamwise 
vortex  rods  lead  to  mixing  enhancement  via  entrainment,  but  also  maintain  high 
product  concentration  near  their  axes  (Fig.  7a). 

In  the  streamwise  sections,  the  product  concentration  chemges  from  a  wavy 
layer  across  the  spem,  into  a  highly  concentrated  distribution  embedded  within 
the  cores  of  the  streamwise  vortices.  The  streamwise  vortex  rods  generate 
mushroom  structures  in  the  braids.  With  the  intensification  of  the  streamwise 
vorticity  in  concentrated  cores.  At  late  stages,  t  -  16.0,  the  field  of  the 
vortices  causes  severe  thinning  of  the  "braids"  joining  neighboring  rods  which 
almost  become  devoid  of  products.  These  transverse  entrainment  currents  are 
best  appreciated  by  simultaneously  examining  Figs.  10a  and  10b,  which  indicate 
that  the  depletion  of  products  occurs  despite  hi^  reaction  rates  in  regions 
separating  the  streamwise  eddies.  Thus,  the  strain  and  entrainment  fields 
tend  to  reshape  the  structure  of  the  product  concentration  such  that  high 
product  concentration  coincides  with  high  concentration  of  vorticity,  and 
falls  rapidly  as  one  moves  into  zones  of  small  vorticity. 

The  effect  of  the  translation  of  the  spanwise  core  on  the  product  and 
reaction  rate  distributions  is  inspected  in  Fig.  9.  At  early  stages,  t  <  8.0, 
combustion  occurs  in  an  almost  uniform  layer  across  the  span.  Initial  growth 
of  the  translative  instability  causes  a  wavy  deformation  of  the  reaction  zone 
and  of  the  region  of  high  product  concentration.  The  distribution  of  the 
latter  is  ccm^licated  by  the  entrainment  of  the  braids,  which  trap  layers  of 
unburnt  fluid  on  the  top  and  bottom  sides  of  the  deformed  core.  At  later 
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stages,  the  .iiaturation  of  the  three-dimensional  instabilities  yields 
distributions  of  even  higher  complexity. 

The  entrainment  of  mushroom  structures  formed  in  the  braids  towards  the 
cores  remains  easy  to  detect.  Within  these  structures,  zones  of  high  product 
coincide  with  the  axes  of  the  corresponding  vortex  rods,  while,  as  expected, 
reaction  rates  remain  high  across  the  sp2m.  The  deformation  associated  with 
the  maturation  of  the  translative  instability  is  harder  to  analyze,  as  the 
core  of  vorticity  to  cross  the  plane  of  the  figure.  The  streamwise  vortices 
induce  a  spinning  motion  which  enhances  the  chemical  reaction  (Fig.  9b)  so 
that  their  aoces  become  zones  of  high  product  concentration.  The  motion 
induced  by  these  vortices,  visualized  by  the  redistribution  of  the  reaction 
zone,  resembles  the  deformation  which  accompanies  the  formation  of  the  vortex 
rods  and  can  be  predicted  by  the  streamwise  vorticity  distribution  shown  in 
Fig.  5.  Interpretation  of  the  product  and  reaction  distributions  shown  in 
Fig.  9  requires  examination  of  several  sections  simultaneously.  Such  an 
exercise  might  be  conducted  by  conparing  the  distributions  of  Figs.  9  and  8, 
with  the  motion  of  material  particles  shown  in  Figs.  2c  and  2b.  The 
comparison  yields  a  clear  illustration  of  the  motion  of  the  core  of  vorticity, 
and  of  the  transverse  diffusion  fluxes  associated  with  the  growth  of  three- 
dimensional  inst2ibilities. 

V.2  REACTING  LAYER  AT  HIGH  DAMKOHLER  NUMBER 

We  now  consider  the  evolution  of  the  reaction  field  for  a  three- 
dimensional  layer  at  high  Oamkohler  and  Peclet  numbers,  Da  -  10,  and  Pe  - 
1000.  As  in  the  previous  section,  the  analysis  is  conducted  by  plotting  the 
product  concentration  and  normalized  reaction  rate  contours  in  the  spanwise 
sections  y  »  3.3  and  1.6  and  the  streamwise  sections  x  -  6.6  and  x  -  2.0,  in 


Figs.  11-14  respectively.  In  each  figure,  snapshots  of  the  correspcndir.g 
variables  are  shown  at  t  =  12.0  and  16.0. 

Comparison  of  figs  Ixa,  12a,  13a  and  14a,  with  their  respective 

counterparts,  Figs.  7a,  8a,  9a,  aind  10a,  reveals  that  the  increase  in  the 

Damkohler  number  does  not  chzmge  the  structure  of  the  product  distribution. 

Despite  high  reaction  rates,  products  are  reorganized  by  the  flow  in  'uch  a 

way  as  to  force  a  correspondence  between  zones  of  high  vorticity  and  high 

product  concentration.  This  similarity  is  in  agreement  with  results  of  two- 
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dimensional  computations. 

At  low  Damkohler  number,  a  small  reaction  rate  confines  the  reaction  zone 
to  the  initial  well-mixed  region.  At  high  Damkohler  number,  the  reaction  zone 
undergoes  a  chemge  as  the  flow  evolves,  from  a  region  inside  the  large 
cores/rods  to  a  thin  zone  around  their  outer  boundaries.  This  is  expected 
since  fast  reaction  rates  lead  to  "conplete"  combustion  in  well-mixed  zones 
and  therefore  to  the  migration  of  the  reaction  zone  to  reactants-rich  regions. 

As  shown  in  fig.  11b,  the  product  concentration  reaches  values  close  to 
one  inside  the  core  of  the  spanwise  eddy.  Similarly,  the  section  located  at  y 
■  1.6,  shown  in  Fig.  12b,  indicates  that  the  reaction  rate  drops  first  at  the 
center  of  the  eddy.  At  late  stages,  the  reaction  zone  migrates  towards  the 
outer  edges  of  the  eddy  core.  Burning  takes  place  in  thin  regions  located 
around  the  outer  edges  of  the  core  of  vorticity^  while  little  or  no  reaction 
occurs  within  the  core.  A  similar  redistribution  of  the  reaction  zone  is 
observed  in  the  spanwise  plane  y  ■  3.3.  However,  this  effect  is  less 
pronoxinced  near  the  spanwise  core,  vdiile  the  braids  continue  to  support  the 
chemical  reaction.  The  difference  between  the  two  sections  is  related  to  the 
presence  of  streamwise  vortices  x4iose  stirring  action  results  in  enhanced 
transverse  diffusion  fluxes. 


The  mechanisms  by  which  the  streamwise  vortices  enhance  the  rate  of 
burning  differ  between  the  core  and  braid  regions.  Within  the  core, 
streamwise  vortices,  generated  by  the  deformation  of  the  core  itself,  cause  a 
spinning  motion  that  is  weakly  affected  by  the  underlying  two-dimensional 
flow.  This  results  in  a  limited  enhancement  in  burning  rates  within  a  large, 
product-dominated  zone.  On  the  other  hamd,  despite  the  presence  of  streamwise 
vortex  rods  within  the  braids,  the  strain  field  associated  with  the  underlying 
two-dimensional  flow  prevents  the  formation  of  a  thick  product  region.  Hence, 
high  product  concentrations  and  high  reaction  rates  coexist  in  this  region. 

The  streamwise  sections  of  Figs.  13b  and  14b  provide  ^ulother  view  of  the 
evolution  of  the  reaction  zone  in  the  spanwise  direction.  At  t  -  12.0,  the 
deformation  of  the  product  interfaces  is  dictated  by  the  growth  of  three- 
dimensional  instabilities.  Combustion  occurs  in  two  thin  layers  enclosing  the 
spanwise  vorticity  core,  and  in  the  braids  which  start  to  roll  under  the 
action  of  the  streamwise  vortex  rods.  At  later  stages,  the  reaction  becomes 
confined  to  the  core  of  the  vortex  rods  and  a  thin  layer  trapped  between  the 
braids  and  the  core  region.  '  The  region  of  lowest  reaction  rate,  located 
within  the  spanwise  core,  corresponds  to  a  zone  of  small  streamwise  vorticity 
separating  neighboring  vortices.  On  the  other  hand,  cross-sections  of  the 
braids  shown  in  Fig.  14b  confirm  our  earlier  claim  regarding  the  competing 
influence  of  the  streamwise  vortex  rods  and  of  the  two-dimensional  strain 
field.  Despite  the  extreme  thinning  of  the  braids  by  the  strain  field  of  the 
spanwise  core,  and  the  continuous  entrainment  of  the  products  of  reaction,  the 
braids  support  high  reaction  rates  and  high  product  concentration  near  the 
axes  of  the  vortex  rods. 


V.3  TOTAL  MASS  OF  PRODUCTS 


The  evolution  of  the  total  mass  of  products  formed,  ri( t)  =  J  o  Cp  dx  =  Z 
D  Cp  ^  dV^,  is  shown  in  Fig.  15  for  all  21  cases  considered.  Figure  15 
contains  three  plots,  corresponding  to  (a)  Pe  -  1000,  (b)  Pe  -  500,  and  (c)  Pe 
=  250,  each  showing  changes  in  M(t)  for  seven  values  of  the  Damkohler  number. 
Da  =  0.1,  0,2,  0.4,  1.0,  2.5,  5.0,  and  10.  Variations  in  the  Peclet  number 
are  found  to  induce  minor  changes  in  the  product  and  reaction  zone  structures. 
Thus,  at  high  Peclet  numbers,  the  latter  are  mainly  governed  by  the  convective 
flow  field.  However,  the  reaction  rates  are  strongly  dependent  on  the  amount 
of  mixing,  and  hence  on  the  coefficients  of  mass  diffusion. 

At  high  values  of  the  Damkohler  number,  the  chemical  reaction  leads  to  eui 

almost  immediate  and  complete  burning  of  the  initial  mixed  region.  Curves 

corresponding  to  Da  >  2.5  are  characterized  by  a  short  and  sharp  initial 

growth  period,  which  does  not  depend  on  the  Peclet  number.  This  is  expected 

since,  for  fast  chemistry,  regions  of  mixed  reactants  bum  inmediately.  An 

abrupt  transition  to  a  controlled  growth  regime  follows,  indicating  that  a 
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diffusion-limited  reaction  is  reached.  The  burning  rates  in  this  regime 

increase  with  decreasing  Peclet  numbers,  as  indicated  by  the  slopes  of  the 
corresponding  M(t)  curves  at  later  stages  between  Figs  15a,  15b  and  15c. 

At  low  values  of  the  Damkohler  number.  Da  <  0.4,  the  chemical  reaction  is 
too  slow  for  such  a  transition  to  occur.  The  evolution  of  M(t)  exhibits  a 
monotonic  increase  throughout  the  duration  of  the  computations.  As  previously 
mentioned,  the  reaction  remains  confined  to  the  initial  mixed  region.  This 
observation  is  supported  by  noting  that  the  mass  of  products  formed  at  the  end 
of  the  computations  at  low  Damkohler  number  is  less  than  that  reached  at  the 
time  of  the  transition  to  diffusion-limited  regime  in  the  high  Damkohler 
number  computations.  This  explains  the  fact  that,  at  late  stages,  the  burning 
rates  at  low  Damkohler  nxjmber  are  higher  than  those  at  high  Damkohler  number. 
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and  that  gains  in  the  amount  of  products  formed  achieved  by  decreasing  the 
Peclet  number  are  smaller  when  the  Damkohler  number  is  small. 


VI  CONCLUSIONS 


In  this  work,  the  evolution  of  reacting  shear  layers  are  numerically 
determined  in  the  high  Peclet  number  regime,  for  a  wide  range  of  Damkohler 
numbers.  Computed  results  show  that  the  product  distribution  is  shaped  by  the 
convective  field  induced  by  spanwise  and  streamwise  vortex  structures  which 
form  due  to  the  growth  of  essential  instability  of  the  flow.  The  rollup  of 
spanwise  vorticity  leads  to  the  creation  of  concentrated  vortex  cores,  and 
braids  which  join  neighboring  cores.  Entrainment  currents  associated  with 
these  structures  force  the  migration  of  products  from  the  braids  towards  the 
cores,  while  their  induced  strain  field  causes  a  severe  thinning  of  the  braids 
and  of  the  reaction  zone  supported  therein.  However,  streamwise  vortices, 
which  are  generated  as  a  result  of  growth  of  three-dimensional  instability  and 
ace  intensified  by  stretch,  significantly  affect  the  flow  at  later  stages, 
resulting  in  substantial  deviation  from  the  two-dimensional  situation. 

The  maturation  of  the  streamwise  vortices  into  strong  streamwise  rods, 
and  the  aitplification  of  the  translative  instability  are  accompanied  by 
spanwise  variations  in  the  reacting  field  and  the  formation  of  mushroom 
stnictures.  Mixing  and  burning  enhancement  is  achieved  through  the  transverse 
entrainment  fluxes.  The  entrainment  fluxes  cause  a  reorganization  of  the 
product  distribution  such  that  zones  of  high  product  concentration  always 
correspond  to  zones  of  hi^  magnitude  of  vorticity.  Thus,  products  tend  to 
migrate  towards  the  core  of  the  spanwise  vortices  and  towards  the  axes  of  the 
streamwise  vortex  rods. 
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Vv'hile  the  product  distribution  is  dictated  by  the  flow,  and  is 
insensitive  to  variation  of  the  Damkohler  number,  the  structure  of  the 
reaction  cone  depends  strongly  on  the  latter.  At  low  Damkohler  number, 
combustion  occurs  in  distributed  zones  located  within  the  cores  of  the 
vortices,  and  is  confined  to  the  initial,  well-mixed  region.  As  the  Deimkohler 
number  increases,  complete  combustion  is  achieved  within  the  cores  of  the 
vortices,  thus  causing  a  migration  of  the  reaction  zone  towards  their  outer 
edges.  The  motion  of  the  reaction  zone  towards  regions  of  higher  strain  rates 
result  in  a  substantial  chauige  in  its  structure,  as  considerable  thinning  of 
the  latter  is  observed.  Extension  of  the  computations  to  study  pairing  among 
several  eddies  and  to  accommodate  high  heat  release,  compressible  flow  models 
with  complex  chemical  reactions,  is  currently  contenplated. 
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FIGURE  OiPTIONS 


Figure  1.  Three-dimensional  perspective  view  of  the  material  surface  initially 
lying  in  the  plane  z  -  0,  and  carrying  the  highest  spanwise  vorticity.  The 
plots  are  generated  at  t  »  12.0  and  16.0  from  the  point  of  view  of  an  observer 
located  at  (48,24,48).  x-  denotes  the  streamwise  direction,  y-  the  spanwise 
direction,  and  z-  is  cross-stream  direction.  A  constant-y  plauie  is  called  a 
spanwise  section,  v^ile  a  constant-x  is  called  a  streamwise  section. 

Figure  2.  Intersection  of  the  Lagrangian  mesh  at  t  -18.0  with  the  planes 
defined  by  (a)  y  -3.3,  (b)  y  -  1.6,  (c)  x  «  6.6  and  (d)  x  -  2.0.  Results  are 
obtained  using  rectangular  transport  elements. 

Figure  3.  Contours  of  the  spaurwise  vorticity,  shown  in  the  x-z  plane 
located  at  y  -  3.3.  The  plots  are  generated  at  t  »  12.0  and  16.0.  Dashed 
lines  are  used  to  represent  negative  values. 

Figure  4.  Contours  of  the  spanwise  vorticity,  shown  in  the  x-z  plane 
located  at  y  -  1.6  at  t  «  12.0  amd  16.0. 

Figure  5.  Contours  of  the  streamwise  vorticity,  shown  in  the  y-z  plauie 
located  at  x  -  6.6  at  t  -  12.0  and  16.0. 

Figure  6.  Contours  of  the  streamwise  vorticity,  shown  in  the  y-z  plane 
located  at  x  -  2.0  at  t  -  12.0  and  16.0. 

Figure  7.  Contours  of  (a)  product  concentration  and  (b)  reaction  rate 
generated  in  the  x-z  plane  located  at  y  -  3.3.  The  plots  are  generated  at  t  - 
ir  0  and  16.0  for  a  reacting  layer  with  Da  ■  0.1  and  Pe  ■  1000. 

Figure  8.  Contours  of  (a)  product  concentration  and  (b)  reaction  rate 
generated  in  the  x-z  plane  located  at  y  •  1.6.  The  plots  are  generated  at  t  - 
12.0  and  16.0  for  a  reacting  layer  with  Da  ■  0.1  and  Pe  ■  1000. 
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Figure  9.  Contours  of  (a)  product 

generated  in  the  y-z  plane  located  at  x 
12.0  and  16.0  for  a  reacting  layer  with 
Figure  10.  Contours  of  (a)  product 

generated  in  the  y-z  plauie  located  at  x 
12.0  and  16.0  for  a  reacting  layer  with 
Figure  11.  Contours  of  (a)  product 

generated  in  the  x-z  plzuie  located  at  y 
12.0  and  16.0  for  a  reacting  layer  with 
Figure  12.  Contours  of  (a)  product 

generated  in  the  x-z  plane  located  at  y 
12.0  and  16.0  for  a  reacting  layer  with 
Figure  13.  Contours  of  (a)  product 

generated  in  the  y-z  plane  located  at  x 
12.0  and  16.0  for  a  reacting  layer  with 
Figure  14.  Contours  of  (a)  product 

generated  in  the  y-z  plane  located  at  x 
12.0  and  16.0  for  a  reacting  layer  with 
Figure  15.  Evolution  of  the  total  mass 
layers  at  (a)  Pe  ■  1000,  (b)  Pe  -  500, 
is  generated  for  Damkohler  nunibers  Da  * 


concentration  and  (b)  reaction  rate 
»  6.6.  The  plots  are  generated  at  t  = 
Da  »  0.1  and  Pe  -  1000. 
concentration  and  (b)  reaction  rate 
*  2.0.  The  plots  are  generated  at  t  - 
Da  -  0.1  and  Pe  -  1000. 
concentration  and  (b)  reaction  rate 
-  3.3.  The  plots  are  generated  at  t  - 
Da  ■  10  and  Pe  -  1000. 
concentration  and  (b)  reaction  rate 
«  1.6.  The  plots  are  generated  at  t  - 
Da  -  10  and  Pe  -  1000. 
concentration  euid  (b)  reaction  rate 
-6.6.  The  plots  are  generated  at  t  - 
Da  ■  10  and  Pe  -  1000. 
concentration  and  (b)  reaction  rate 
-2.0.  The  plots  are  generated  at  t  - 
Da  »  10  and  Pe  »  1000. 
of  products,  M(t),  for  reacting  shear 
and  (c)  Pe  -  250.  In  each  plot,  M(t) 
0.1,  0.2,  0.4,  1,  2.5,  5,  and  10. 
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ABSTRACT 


1.  INTRODUCTION 


The  transport  element  method,  a  Lagrangian,  grid-free 
scheme  based  on  the  extension  of  the  vortex  method  to 
reacting  compressible  flow  at  low  Mach  number,  has  been 
applied  to  simulate  a  spatially  developing,  reacting  shear  layer 
with  unmixed  reactants.  The  method  employs  computational 
elements  which  move  along  particle  trajectories  and  transport 
vorticity  and  scalar  gradients.  In  the  case  of  a  reacting  flow, 
vorticity  changes  due  to  density  variation,  scalar  gradients 
change  due  to  the  chemical  reaction,  and  volumetric  expansion 
adds  an  expansion  field.  Solutions  are  obtained  for  a  forced 
shear  layer  at  different  Damkohler  numbers  and  enthalpy  of 
reaction  to  study  the  effect  of  combustion  heat  release  rate  on 
the  development  of  the  large  scale  structures.  Forcing  is  used 
to  ensure  roll-up  within  the  computational  domain.  We  show 
that  heat  release  enlarges  the  size  of  the  fundamental  eddies, 
stretching  their  streamwise  dimension  and  slightly  reducing 
their  cross  stream  dimension,  while  their  overall  size  remains 
almost  the  same.  Along  with  forcing  at  the  fundamental  and 
subharmonic  frequencies,  heat  release  increases/decreases  the 
size  of  the  in-phase/out-of-phase  eddies.  The  non-uniform 
acceleration  of  the  eddies  in  the  streamwise  direction  causes 
their  relative  locations  to  deviate  from  that  of  a  uniform-density 
layer  and  thus  modifies  the  pairing  process  into  a 
tearing/gulping  process.  Results  also  show  that  the  enthalpy 
of  reaction  is  more  important  than  the  reaction  frequency  factor 
in  affecting  the  flow  dynamics.  For  the  same  parameters,  the 
variable-density  layer  grows  slower  than  its  uniform-density 
counterpart. 
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Reacting  shear  flow  at  high  Reynolds  number  is  of 
fundamental  importance  to  the  physics  of  combustion  and  to  a 
wide  range  of  engineering  applications.  Due  to  the  complexity 
of  the  problem,  numerical  simulation  has  been  recognized  as 
an  important  tool  in  studying  the  different  phenomena 
governing  turbulent  reacting  flows  [1,2].  To  perform  these 
simulations  accurately,  we  have  developed  a  Lagrangian- 
based,  grid-free,  adaptive  method  to  simulate  reacting  shear 
flow  and  applied  this  methodology  to  several  configurations  of 
fundamental  and  practical  significance.  This  method,  the 
transport  element  method,  has  been  developed  by  extending 
the  vortex  element  method  so  that  compressible  flow  effects 
and  scalar  transport  can  be  included  in  the  simulations.  So  far, 
the  focus  of  the  application  of  this  method,  and  most  other 
numerical  methods,  has  been  the  idealized,  temporally 
evolving  layer  in  which  one  assumes  that  a  reference  frame 
moving  at  the  phase  speed  of  the  fundamental  structures 
follows  the  evolution  of  the  layer  [3-6].  This  idealization, 
however,  leaves  out  several  important  effects,  which  the 
spatially  developing  shear  layer  possesses,  that  have  been 
detected  experimentally  and  recognized  in  modeling  analysis. 
The  purpose  of  this  paper  is  to  apply  the  transport  element 
method  to  a  spatially  developing  shear  layer  and  to  investigate 
the  effect  of  heat  release  on  the  fluid  flow  in  this  problem. 

Spatially  developing  shear  layer  exhibits  an  asymmetric 
entrainment  behavior  which  tends  to  drift  the  composition  of 
the  large  scale  structures  towards  the  species  convected  by  the 
high-velocity  stream  [7,8].  If  the  two  streams  have  different 
densities,  or  different  momenta,  that  also  contributes  to  the 
exact  composition  of  the  mixed  fluid  inside  the  large  eddies  in 
ways  similar  to  that  of  the  velocity  ratio  although  through  a 
different  mechanism  [9].  We  found  in  previous  simulations 
that  forcing  can  alter  this  composition  in  a  way  which  depends 
on  the  momentum  ratio  between  the  two  streams.  Most  of 
these  findings  were  supported  by  linear  stability  analysis  and 
several  experimental  studies.  While  is  is  possible  to  capture 
the  effect  of  the  density  ratio  on  the  change  of  Ac  composition 
of  the  large  structures  in  a  temporal  simulation,  it  is  not 
possible  to  study  the  effect  of  the  velocity  ratio,  or  the 
combined  effects  of  Ae  velocity  and  density  ratios,  in  such  a 
configuration.  This  is  because  a  temporal  simulation  assumes 
perfect  dynamic  symmetry  across  the  mean  velocity  contour 
and  any  sub^uent  changes  in  Ae  layer  structures  are  limited 
by  Ais  restriction.  Our  analysis  of  the  spatially  developing 
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variable-density  shear  layer  showed  that  the  results  of  the 
temporal  analysis  are  either  limited  to  a  narrow  range  of 
parameters  or  are  literally  wrong  if  the  effect  of  forcing  is 
considered.  We  conclude  that  it  is  necessary  to  perform 
spatially  developing  shear  layer  simulations  if  the  physics  of 
the  reacting  flow  are  to  be  properly  modeled. 
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+  pVu  =  0 
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The  simulations  of  reacting,  spatially  developing 
mixing  layers  with  uniform  density,  i.e.,  when  the  dynamic 
effects  of  heat  release  were  neglected,  have  been  used  before 
to  discern  the  structure  of  the  reaction  zone  and  the  effects  of 
the  Damkohler  number  on  the  relative  location  of  the  reaction 
zone  with  respect  to  the  large  structures.  It  was  concluded  that 
at  small  DatrJcohler  number,  the  reaction  is  most  intense  near 
the  center  of  the  large  eddies,  while  as  the  Damkohler  number 
increases,  the  reaction  zone  moves  outwards  towards  the  outer 
edges  of  the  eddies.  The  exact  thickness  of  the  reaction  zone, 
which  resembles  that  of  a  strained  diffusion  flame,  is  a 
function  of  the  Damkohler  and  Reynolds  numbers.  It  was  also 
found  that,  under  conditions  of  unity  stoichiometry,  a  strong 
similarity  exists  between  the  products  concentration  field  and 
the  vorticity  field,  suggesting  a  possible  modeling  approach 
for  this  problem.  In  all  these  simulations,  the  flow  was 
restricted  to  uniform  density.  The  results  of  the  spatially 
developing,  nonreacting,  variable  density  shear  layer  showed 
that  even  without  volumetric  expansion,  density  variation  can 
influence  the  layer  dynamics  in  a  nontrivial  way.  Thus,  it  is 
important  to  introduce  the  dynamic  effects  of  heat  release  into 
the  spatially  developing  layer  simulations. 

Dynamic  effects  of  heat  release  have  been 
demonstrated  before  in  a  temporally  developing  shear  layer. 
There  are  two  controversial  issues  regarding  the  effect  of 
combustion  on  the  flow  dynamics;  (a)  whether  heat  release 
reduces  the  growth  rate  of  the  layer  via  baroclinic  vorticity 
generation  or  via  volumetric  expansion;  and  (b)  whether  there 
is  actually  a  substantial  reduction  in  the  size  of  the  eddies  or  it 
is  just  a  delay  in  the  instability  development.  Although  it  has 
been  established  that  heat  release  at  the  early  stages  can  reduce 
the  growth  rate  of  the  layer,  it  was  also  found  that  forcing 
beyond  the  linear  range  can  be  used  to  overcome  this  effect. 
This  was  demonstrated  in  the  temporal  layer  only  and  will  be 
investigated  in  the  spatial  case  in  this  paper.  Deciding  on  the 
separate  role  of  baroclinic  vorticity  and  volumetric  expansion 
requires  simulations  in  which  one  mechanism  is  disabl^  while 
the  other  is  active.  This  will  be  the  subject  of  future  work. 
The  second  issue  could  not  be  resolv^  by  the  temporal 
simulations  only  since  the  behavicn  of  the  spatially  developing 
layer  differs  from  that  of  the  temporal  layer  especially  in  the 
pairing  mode  which  dominate  the  activities  in  the  downstream 
section. 


In  this  paper,  we  present  simulations  of  a  spatially 
developing  shear  layer  and  use  our  results  to  answer  some  of 
these  questions.  We  introduce  the  extension  of  the  transport 
element  method  to  a  reacting  compressible  flow  and  its 
application  to  the  spatially  growing  shear  layer.  We  then 
present  results  of  a  nonreacting  and  reacting  shear  layer  at 
different  rates  of  heat  release. 


n.  FORMULATION  AND  NUMERICAL  SCHEME 
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where  V  =  is  the  gradient  operator,  ^  ^ 


the  material  derivative,  p  is  the  mixture  density,  t  is  time, 

M  =  (m,v)  is  the  velocity  vector  in  an  (x,y)  right-handed 
coordinate  system  and  p  is  pressure.  The  m^el  requires  the 
transport  of  three  scalars  sj,j  =1^,3  :  the  temperature,  T,  and 
the  reactants'  mass  fractions,  Y]  and  Y2.  The  product  mass 

fraction  is  Tp  =  1  -  (Ti  +  T2).  The  reaction  rate  is 


T 

w  =  AfpYi  Y2  exp{-  -^)  where  Af  is  the  non-dimensional 

frequency  factor  and  Tq  is  the  normalized  activation  energy. 
For  temperature  tfansport  Pej  =  Pe,  the  Peclet  number,  and  Qj 

=  20  where  Q  is  the  enthalpy  of  reaction  and  0  is  the  mass 
stoichiomerty.  For  the  reactants  Pej  =  PeLe  where  Le  is  the 

Lewis  number.  For  T;,  Qj  =-0  and  for  Y2,  Qj=-1-  Both  the 
Lewis  number  and  the  stoichiometry  are  assumed  equal  to 
unity. 


The  vorticity  equation  is  obtained  by  taking  the  curl  of 
equation  (2); 

"  (5) 

and  the  scalar  gradient  equation  by  taking  the  gradient  of 
equation  (3): 


where  o>  = - is  the  vorticity,  gj  =  'Vsj  the  scalar 

dx  dy 

gradient  and  /I  is  a  unit  vector  normal  to  the  plane  of  motion. 
The  density  gradient  is  obtained  from  the  temperature  gradient 
using  equation  (4)  and  the  pressure  gradient  from  the  velocity 
field  using  equation  (2).  *1110  velocity  is  decomposed  into  a 
vorticity  induced  solenoidal  component  and  two  iirotational 
components;  one  induced  by  volumetric  expansion  and  the 
other  by  the  boundary  conditions  (u.n  prescribed  at  the 
boundaries,  n  denoting  the  unit  vector  normal  to  the 
boundary),  i.e.. 


We  consider  the  motion  of  a  two-dimensional, 
compressible,  chemically  reacting  flow  at  low-Mach  number 
and  high  Reynolds  number.  We  assume  that  both  reacunts 
and  products  behave  as  perfect  gases  with  equal  molecular 
weight  and  constant  transport  properties,  that  the  reaction  is 
sin^e  step  and  the  kinetics  are  of  the  Arrhenius  type.  Under 
these  assumptions  the  governing  equations  in  non-mmensionaJ 
form  are: 
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where  Ua,  =  Vx(v*),Mi  =  V0e,M^=V0>,  andv'.^e,^ 
satisfy; 
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and  in  the  second  step  we  add  the  contribution  of  the  baroclinic 
torque  while  approximating  the  material  acceleration  of  the 
particle  by  a  two-step  iteration  forward  difference  scheme. 


The  numerical  integration  of  the  vorticity  transport 
equation  starts  by  writing  the  vorticity  as  a  summation  over  the 
fields  of  a  finite  number  of  vortex  elements  of  finite 
overlapping  circular  cores: 
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Note  that  the  evolution  of  the  vorticity  field  requires 
knowledge  of  the  density  field  which  is  coupl^  to  the  former 
via  the  mechanisms  of  bMoclinicity  and  volumetric  expansion. 
The  density  field  is  obtained  from  the  tempierature  field  via  a 
discrete  equivalent  of  equation  (4).  Similarly,  the  density 
gradient  is  obtained  by  taking  the  gradient  of  equation  (4). 
The  expansion  velocity  Ue  can  be  obtained  using  equation  (9) 
in  a  similar  fashion  to  the  calculation  of  the  vortical  velocity 
from  the  vorticity  field,  i.e.. 


where  Hit)  =  6^(0  is  the  circulation  of  the  ith  element,  ^ 

^iU),  XiU)y  and  hfu)  denote  its  vorticity,  location,  and  area,  =  X  p  ^G^x-Xi(t)) 

respectively,  and 
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fdr)  = 

JtS^  (12) 

is  the  second-order  Gaussian  core  function  with  core  radius  S 
used  to  obtain  a  second-order  discretization  of  the  vorticity 
field 


where 

VGs(I)  =  ^  a-exp(-  4)) .  »■  =  ® 

2;ir2  ^2 

is  the  desingularized  gradient  of  the  kernel  of  the  Poisson 
equation. 


The  vortical  velocity,  Uta,  induced  by  the  vorticity  field 
in  equation  (1 1)  is  given  in  terms  of  the  desingularized  Biot- 
Savart  law: 


ujZt)  =  X 
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where 


The  integration  of  the  scalar  gradients  transport 
equation  (6)  is  initiated  by  discretizing  the  gradient  fields 
amongst  a  finite  number  of  transport  elements  in  a  similar 
fashion  to  the  discretization  of  the  vorticity  field,  i.e., 

y  - 
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where  C/i(0  =  l;i(0  *?(f),  is  the  weighted  scalar  gradient  of  the 
ith  element  and  gji  the  corresponding  gradient. 


is  the  desingularized  kernel  of  the  Biot-Savart  integral  and 
r  =  etl. 

The  evolution  of  the  vorticity  field  is  found  by 
updating  the  locations  of  the  vortex  elements,  and  their 
associated  vorticity  values  according  to  equation  (5).  The 
locations  of  the  elements  are  determined  by  advecting  the 
vortex  elements  with  the  local  velocity  vector  along  particle 
trajectories.  This  is  done  by  integrating: 


The  scalar  fields  can  be  reconstructed  from  their 
gradients  via  a  convolution  over  the  discretized  field: 

N  -  _ 

Sjixj)  =  X  C><(')  •  ^Gi0-Xi{t))  +  Spj(x,t) 

<=i  (21) 

where  Spj{x,t)  is  a  component  added  to  satisfy  the  boundary 
conditions. 


with  initial  condition  XiiZifi)  =  i.e.  Xi  ihc  Lagrangian 

coordinate  of  the  particle  originally  located  at  X  j.  Meanwhile, 
the  evolution  of  the  vorticity  associated  with  the  elements  is 
established  by  considering  the  latter's  circulation  which 
evolves  according  to: 

Dfj(f)_  Vpi^Dui 

Of  Pi  Of  (jgj 

Equation  (16)  is  integrated  in  two  fractional  steps.  In  the  first 
$t^,  the  convection  step,  we  have: 


To  simulate  the  evolution  of  the  scalar  gradient  fields, 
the  transport  elements  are  advected  in  a  similar  fashion  to  the 
vonex  elements  and  their  gradients  are  updated  according  to 
equation  (6)  which  is  integrated  in  two  fractional  steps.  In  the 
first  step,  the  convection  and  reaction  step,  we  consider: 


Rather  then  integrating  numerically  equation  (22)  an  alternative 
approach  is  considered.  The  scalar  g^ients  are  related  to  the 
Lagrangian  deformation  of  elementary  surface  elements  whose 

normal  coincides  with  the  direction  of  the  gradient.  If  S 
denotes  such  an  elementary  surface,  then  the  evolution  of  its 
absolute  value  is  governed  by: 
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dt 


(23) 


where  n  is  a  unit  vector  lying  in  the  direction  of  the  scalar 
gradient.  By  considering  the  dot  product  of  n  with  equation 
(22)  and  utilizing  equation  (23)  it  can  be  shown  that: 

dI  UjI  l^>l  \ 
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Equation  (3)  also  suggests  that 


B^  =  Qjd^Ssj 
Dt  ‘  ’ 


ds, 
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where  dsj  is  the  variation  of  the  scalar.  Comparison  of 
equations  (24)  and  (25)  implies  that: 
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Assuming  negligible  directional  effects  of  diffusion  on  the 
scalar  gradients  the  constant  C  depends  only  on  initial 
conditions.  Thus,  knowledge  of  this  constant,  of  the  material 
surfaces  (which  are  readily  available  in  this  Langrangian 

scheme  and  should  yield  both  n  and  S  )  and  of  Ssj  (via 
numerical  integration  of  equation  (25))  establishes  the  time 
evolution  of  the  scalar  gradient  due  to  convection  and 
reaction. 


In  the  second  fractional  step  we  add  the  effect  of 
diffusion: 
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This  is  done  by  expanding  the  cores  of  the  elements  according 
to: 


Pet 


(28) 


where  So  is  the  ewe  radius  at  time  zero. 

Previous  application  of  the  vortex  and  transport 
element  method  to  the  spatially-developing,  uniform  density 
shear  layer  has  been  used  to  validate  the  scheme  by  comparing 
the  numerical  results  for  velocity  and  scalar  concentration 
statistics  to  experimental  measurements. 


ra.  FLOW  GEOMETRY  AND  BOUNDARY  CONDITIONS 


The  reacting  shear  layer  is  formed  by  the  merging, 
downstream  of  a  semi-infinite  splitter  plate  of  vanishingly 
small  thickness,  of  two  streams  of  reactant  fluids  of  unequal 
velocity  (x=0).  The  two  fluids  which  are  assumed  perfect  and 
with  the  same  tixilecular  weight  are  at  the  same  pressure  and 
temperature,  thus,  same  density.  The  top  stream  has  velocity 
Ui  and  consists  of  reactant  Yj,  while  the  bottom  stream  has 
velocity  U2  and  consists  of  reactant  Yj.  Variables  are 
normalized  with  respect  to  Ui,  Tg,  and  H,  the  top-stream 


velocity,  the  ambient  temperature,  and  the  channel  height, 
respectively.  The  solution  domain  consists  of  the  channel 
region  lying  downstream  of  the  splitter  plate  and  bounded  at 
^max,  0  <x  <  Xmax.  0  <  y  <  H.  The  top  and  bottom  walls 
are  modelled  as  rigid,  slip,  impermeable  planes,  so  that  the 


<7J; 

boundary  conditions  v  =  0  and  =  0  at  y  =  0  and  y  =  H  are 

imposed,  where  sj  is  any  of  the  three  transported  scalars  T, 
Yl.  ¥2-  These  conditions  are  satisfied  by  conformally 
mapping  the  entire  channel  region  on  the  upper  half  of  a 
complex  plane,  and  using  the  appropriate  image  system  of  the 
transport  elements.  At  the  downstream  section,  a  condition  of 
vanishing  vorticity  and  scalar  gradient  is  used  as  outflow 
boundary  condition,  and  is  applied  by  removing  the  transport 
elements  which  cross  the  x  =  Xmax  plane.  At  the  inlet  section, 
the  vorticity,  velocity,  top  reactant  mass-fraction-gradient  and 
mass-fraction  prof’Ies  are  specified  as  follows: 
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where  erf  denotes  the  error  function.  The  bottom  reactant's 
mass-fraction-gradient  and  mass-fraction  profiles  are  mirror 
images  of  the  top  reactant’s  profiles.  In  equations  (29-32),  o 
denotes  the  standard  deviation  of  the  Gaussian  profile  which 
describes  the  initial  vorticity  profile.  Thus  the  mass-fraction 
profiles  are  thinner  than  the  vorticity-velocity  profiles  and  are 
displaced  about  the  centerline  of  the  channel  by  a  distance  dy. 

The  inlet  product  profiles  can  be  obtained  from  the  reactant 
profiles  via  Yp  =J-(Yj  +  Y2)-  The  temperature  and  product 
profiles  are  assumed  to  be  similar  thus  the  latter  are  used  to 
define  the  former. 


— (x=0,y,f)  =  Tf,-^x=0,y,t) 
dy  dy 

T{x=0,y,t)  +  Tfi  Ypix=0,y,t) 

where  Tjj  is  the  adiabatic  flame  temperature.  The  thickness  of 
the  vorticity  layer,  2  a,  is  scaled  so  that  the  height  of  the 
channel,  H,  equ^s  twice  the  wavelength  of  the  most  unstable 
fundamental  mode  of  the  initial  vorticity  distribution  of  the 
uniform-density  layer.  For  the  profiles  considered  here,  2ar  = 
0.08  H. 

Initially  the  layer  is  assumed  to  be  flat  and  the  above 
inlet  conditions  are  assumed  valid  throughout  the  domain.  The 
vorticity  and  scalar-gradient  fields  are  discretized  by 
distributing  transport  elements  over  nine  material  surfaces 
(lines)  lying  within  the  region  of  significant  vonicity  and 
gradient  variations.  The  layers  are  equally  separated  with  a 
separation  distance  dy  =  0.0234,  yielding  a  total  thickness  of 

the  discretization  domain  4  =  0.1875H  (i.e.  ^  >  2a).  The 
initial  streamwise  separation  distance  between  neighboring 
elements  lying  on  the  same  material  line,  dx  =  dy,  so  that  a 
square  Lagrangian  mesh  is  obtained  (Fig.  I.).  The  vorticity 
and  scalar-gradient  values  associated  with  the  transport 
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elements  are  found  by  matching  the  discretized  and  assumed 
profiles  at  their  centers,  i.e.,  by  solving  the  linear  system 
obtained  by  equating  the  left-hand  side  of  equation  (1 1)  to  that 
of  equation  (29)  (or  that  of  equation  (20)  to  (31)  or  (33)), 

evaluated  at  X The  value  of  the  core  radius  is  found  by 
iteration  over  those  values  which  ensure  strong  initial  overlap 
so  as  to  minimize  the  integral  error  between  the  assumed  and 
discretized  profile.  For  the  assumed  vorticity  profile  and  for 
the  specified  mesh,  this  procedure  gives  S  -  0.0273.  This 
collocation  scheme  was  found  to  yield  more  accurate  results  on 
the  initial  growth  of  the  instability  waves. 

The  severe  stretching  of  the  La^angian  mesh  used  to 
discretize  the  vorticity  and  scalar  gradients,  which  increases 
the  distance  between  neighboring  elements,  may  lead  to  the 
deterioration  of  the  accuracy  of  the  discretization.  To 
overcome  this  problem,  a  scheme  of  local  mesh  refinement  is 
adopted  whereby  transport  elements  are  continuously 
introduced  and  deleted  to  ensure  overlap  of  the  elements.  In 
our  computation,  strong  overlap  is  enforced  near  the  inlet  of 
the  domain  by  allowing  a  small  maximum  separation  between 
neighboring  elements.  This  condition  is  relaxed  as  we 
approach  the  exit  of  the  domain  by  increasing  this  threshold 
value,  thus  allowing  for  efficient  computations  without 
compromising  the  accuracy  of  the  numerical  scheme. 


IV.  RESULTS 


The  model  described  in  Section  III  has  been  used  to 
study  the  effect  of  heat  release  on  the  dynamics  of  the  spatially 
growing  shear  layer  and  the  influence  of  that  on  the  rate  of 
product  formation.  In  previous  studies  of  temporally  growing 
shear  layers  and  jets,  it  has  been  found,  especially  in  an 
nonpremixed  shear  layer,  that  heat  release  reduces  the  growth 
rate  of  the  instability  and  slows  down  the  roll-up  in  the  linear 
range.  It  has  been  postulated  that  the  volumetric  expansion  in 
the  reaction  zone,  which,  in  the  early  stages,  lies  within  the 
vorticity  layer,  leads  to  the  thickening  of  the  vorticity  layer  and 
thus  slows  down  the  roll-up.  It  has  also  been  established  by 
our  previous  studies  on  a  reacting  jet  that  forcing  at  relatively 
large  amplitudes  can  be  used  to  overcome  this  cfTcct.  If  the 
layer  is  forced  into  its  nonlinear  stages  via  large  amplitude 
oscillations,  the  effect  of  heat  release  on  the  growth  rate  of  the 
large  eddies  is  unnoticcable.  Thus,  in  this  work,  we  use  large 
amplitude  forcing  at  the  splitter  plate  to  initiate  the  roll-up 
before  any  measurable  heat  release  and  thus  avoid  the 
suppression  of  the  instability  due  to  exothermic  energy.  The 
focus  of  the  current  study  is,  therefore,  on  the  effect  of  heat 
release  on  the  large-scale  structure. 

The  forcing  function,  shown  in  Fig.  2,  consists  of  a 
fundamental  frequency  and  its  subharmonic  frequency.  The 
former  is  chosen  to  correspond  to  the  most  unstable  mode  of 
the  uniform  density,  spatially  developing  shear  layer.  It 
should  be  noted  that  it  is  not  possible,  within  the  framew(»k  of 
the  linear  stability  thec^,  to  find  the  most  unstable  mode  in  the 
variable-density,  reacting  shear  layer  since  the  profiles  change 
with  time  as  energy  is  deposited.  The  amplitude  of  both 
modes  is  taken  to  be  the  same  and  is  0.015.  Note  that,  as  will 
be  seen,  forcing  at  such  high  amplitude  affects  the  layer 
structure  in  a  fundamental  way.  We  observe  that  due  to  the 
presence  of  two  frequencies,  one  being  the  subharmonic  Of  the 
other,  two  eddies  form  and  then  pair.  The  eddy  which  forms 
during  the  pan  of  the  cycle  in  which  the  fun^mental  and  the 
subharmonic  components  are  in-phase  is  bigger  than  that 
which  forms  in  the  second  pan  of  the  subcycle  during  which 
the  two  components  are  out-of-phase.  The  difference  between 
the  size  of  the  two  eddies  increases  as  the  amplitude  of  forcing 
increases.  As  we  wiU  show,  this  has  an  imponant  effect  of  the 
pairing  dynamics  of  the  reacting  shear  layer. 


We  have  performed  several  runs  in  which  the  velocity 
ratio  was  kept  constant,  r  =  1I2.  The  Datnkohler  number  was 
maintained  low  by  keeping  Aj  =  50-100,  while  the  heat  release 
was  maintained  high  by  using  Q  =  4-6.  The  non-dimensional 
activation  energy  was  kept  at  Ta=10  in  all  runs.  Calculations 
of  the  reacting  layer  were  also  performed  for  a  uniform-density 
flow,  i.e.,  when  heat  release  is  not  allowed  to  affect  the  flow 
dynamically,  to  isolate  the  effect  of  density  variation.  Keeping 
the  reaction  frequency  factor  low  limits  the  reaction  zone 
structure  to  that  of  a  distributed  zone  in  which  most  of  the 
exothermic  energy  is  deposited  inside  the  large  structures.  We 
have  shown  before  that  increasing  the  Damkohler  number 
causes  the  reaction  zone,  which  starts  inside  the  cores  of  the 
large  structures,  to  move  towards  the  outer  edges  of  the  eddies 
and  to  become  a  thin  reaction  zone  of  the  flamelet  type. 
Limiting  the  simulations  to  low  Damkohler  numbers,  which  is 
not  the  general  case,  is  representative  of  the  early  stages  of 
development  of  the  reacting  layer  at  an  arbitrary  Damkohler 
number  since  the  reaction  starts  always  in  the  well-mixed 
zone,  i.e.,  inside  the  cores  of  the  structures.  Thus,  we  believe 
that  the  conclusions  of  the  study  apply  to  a  wider  range  of  Dg. 

Figure  3  shows  a  comparison  between  the  reacting 
shear  layer  without  and  with  density  variation,  in  both  cases 
with  temperature  variation  due  to  the  chemical  reaction.  The 
computations  were  performed  with  Af  =  50  and  Q  =  6.  The 
shear  layer  is  depicted  using  all  the  vortex  elements  and  their 
velocity  vectors  at  three  time  steps.  The  velocity  is  plotted 
with  respect  to  the  mean  velocity  across  the  shear  layer  to 
highlight  the  boundaries  of  the  large  eddies  better.  We  note 
that  density  variation,  if  coupled  to  temperature  variation  via 
the  equation  of  state,  leads  to  the  generation  of  vorticity  via  the 
baroclinic  effect  and  to  an  additional  velocity  field  by  the 
volumetric  expansion.  The  effect  of  volumetric  expansion  is 
clear  in  the  larger  velocity  vectors  in  the  downstream  sections 
of  the  shear  layer.  The  cross  sections  of  the  eddies  in  the 
variable-density  shear  layer  are  larger  than  in  the  uniform- 
density  layer.  If  the  cross  section  of  the  eddy  is  modeled  by 
an  ellipse  with  a  major  axis  in  the  streamwise  direction  and  a 
minor  axis  in  the  cross-stream  direction,  the  figure  shows  that 
heat  release  causes  the  eddies  to  expand  in  the  streamwise 
direction  so  that  their  major  axis  increases.  On  the  other  hand, 
the  minor  axis  of  the  variable  density  eddies  decreases  slightly 
giving  the  appearance  that  the  roll-up  of  the  fundamental 
^dies  has  bran  suppressed. 

The  expansion  of  the  eddies  in  the  streamwi  •• 
direction,  and  the  overall  acceleration  of  the  structures  as  the 
flow  moves  downstream,  brings  the  large  eddies  closer  than  in 
the  case  of  the  uniform-density  case.  This  has  the  effect  of 
changing  the  phase  relationship  anxmg  the  neighboring  eddies 
from  that  imposed  by  the  upstream  forcing,  i.e.,  the  relative 
locations  of  the  large  eddies  are  different  than  that  in  the 
uniform-density  layer.  Since  pairing  dynamics  depend 
strongly  on  the  phase  relationship  among  neighboring  eddies, 
this  change  of  phase  due  to  energy  release,  as  expected,  leads 
to  modifying  the  pairing  dynatitics.  This  is  seen  in  Fig.  3  as 
the  eddies  in  the  unifexm-density  case  pair  and  form  structures 
which  protrude  significantly  in  the  cross-stream  direction 
while  in  the  variable-density  case,  the  structures  formed  by 
pairing  seem  to  remain  flat.  The  mechanisms  by  which  the 
two  eddies  merge  in  the  variable-density  case  seem  to  be  a 
merging  or  "gulping"  more  than  pairing.  It  also  appears  that 
small  e^ies  are  tom  between  their  larger  neightors.  The 
reduction  of  the  cross-stream  growth  of  the  layers  can  be 
attributed  in  part  to  the  larger  phase  velocity  of  the  eddies  due 
to  the  volumetric  expansion. 

Figure  3  shows  that  although  heat  release  contributes 
to  the  growth  of  the  structures,  it  does  so  primarily  in  the 
streamwise  direction  and  thus  the  growth  of  the  layer,  as 
measured  by  the  cross  stream  divergence  of  the  boundaries  of 
the  vorticity  layer,  seem  to  be  leduc^  by  heat  release.  This  is 


5 


confirmed  in  Fig.  4  where  we  compare  the  results  of  three 
cases  for  which  Af  =  50  and  100  atQ  =  4,  and  Af=  50  aiQ  = 
6  In  these  figures  we  distinguish  between  the  effect  of  the  rate 
and  the  amount  of  energy  deposition  due  to  combustion.  In 
the  first  and  last  cases,  we  observe  the  peculiar  behavior  that 
the  difference  between  the  sizes  of  the  two  eddies  which  form 
during  one  cycle  of  subharmonic  forcing  becomes  larger. 
Note  that  during  each  cycle  of  subhaimonic  forcing,  two 
subcycles  of  fundamentals  exist:  one  in-phase  and  one  out-of¬ 
phase  with  the  subharmonic.  Even  in  the  uniform-density 
layer,  the  eddy  forming  during  the  in-phase  subcycle  is  larger 
than  the  eddy  which  forms  during  the  out-of-phase  subcycle. 
With  heat  release,  the  "in-phase"  eddies  become  larger  and  the 
"out-of-phase"  eddy  becomes  smaller.  This  imbalance  is 
clearly  seen  in  Fig.  4  where  in  the  downstream  part  of  the 
domain,  the  small,  out-of-phase  eddies  seem  to  get  gulped  by 
the  large,  in-phase  eddies. 

The  reason  why  heat  release  contributes  to  the  faster 
growth  of  the  in-phase  eddy  and  the  slower  growth  of  the  out- 
of-phase  eddy  is  not  clear  at  this  point  and  requires  further 
investigation.  What  is  clear  is  that  as  the  rate  of  energy 
release,  as  determined  by  Af,  or  the  amount  of  energy 
released,  as  defined  by  Q,  increase,  the  difference  between  the 
size  of  the  two  eddies  becomes  more  pronounced.  In  both 
cases,  near  the  downstream  end  of  the  domain,  the  small 
eddies  disappear  faster  as  they  are  either  tom  between  the  two 
larger  neighboring  eddies  or  gulped  by  either  of  them.  On  the 
other  hand,  as  more  heat  is  deposited  inside  the  large  eddies, 
they  expand  both  in  the  streamwise  and  cross  stream 
directions.  This  expansion  is  exhibited  by  the  faster  growth  of 
the  layer  with  higher  heat  release  than  with  the  lower  heat 
release. 

The  effect  of  heat  release  on  the  physical  size  of  the 
layer  is  depicted  in  Fig.  3  where  we  show  (a)  the  product 

thickness,  defined  as  \y(Yp'  =  0.01)  -y(  Yp*=  0.01  )\,  and  (b) 
the  reaction  thickness,  defined  as  lyfW  *  0.01  w„ax)  -yi^*  - 
0.01  WmaxlU  both  averaged  over  100  time  steps,  where 
superscripts  and  -  indicate  increasing  and  decreasing  values 
of  the  variables,  respectively.  These  plots,  especially  the 
uniform-density  and  low-heat-release  cases,  exhibit  the  trend 

previously  found  in  forced  shear  layers:  the  thickness  first 
increases  due  to  the  formation  of  the  fundamental  eddy,  stays 
flat  where  the  structures  are  essentially  convected  without 
interaction,  then  increases  again  during  pairing.  However,  as 
the  rate  of  heat  release  increases,  the  size  of  the  structures 
continues  to  increase  beyond  the  first  stage  due  to  the 
volumetric  expansion  resulting  from  the  heat  release  inside  the 
structures.  Moreover,  the  sudden  transition  where  pairing 
starts  in  the  uniform-density  case  disappears.  Instead,  the 
thickness  seems  to  keep  its  gradual  rise  into  the  Aird  stage. 
This  corroborates  our  previous  observation  that  pairing  in  the 
uniform  density  case  is  replaced  by  a  gradual  process  of 
tearing/gulping  that  leads  to  the  development  of  the  large 
eddies. 

Plots  of  the  layer  thickness  show  that  as  the  rate  of  heat 
release  increases,  the  size  of  the  layer  grows.  Since  the 
growth  of  the  layer  is,  to  a  large  extent,  due  to  volumetric 
expansion,  one  can  infer  that  as  the  layer  thickness  grows 
faster,  the  amount  of  product  formed  also  increases.  It  is 
interesting  to  note  that  while  small  rates  of  heat  release  can 
reduce  the  size  of  the  structures  and  delay  the  overall  growth 
of  the  layer  with  respect  to  the  uniform  density  layer,  as  the 
rate  of  heat  release  increases  both  the  overall  size  of  the  layer 
and  the  rate  of  product  formation  become  higher.  It  should  be 
emphasized  that  this  is  only  true  when  high  amplitude  forcing 
is  applied.  When  we  did  not  force  the  layer,  we  found  that 
increasing  the  heat  release  rate  suppressed  the  instability 
leading  to  the  roll-up,  and  the  rate  of  product  formation  was 
reduced.  It  is  probable  that  stronger  subharmonic  forcing 


would  have  led  to  pairing.  However,  investigating  this 
possibility  is  left  for  future  work. 

The  detailed  structure  of  the  reacting  layer  is  shown  in 
Fig.  6  for  the  uniform-density  case  and  Figs.  7  and  8  for  the 
variable-density  cases.  In  these  figures,  the  results  within  two 
sections  of  the  computational  domain  are  shown  in  terms  of 
the  product  mass-fraction,  the  reaction  rate,  and  the  vorticity. 
For  comparison.  Fig.  6  shows  the  product  mass-fraction, 
reaction  rate  and  vorticity  for  the  uniform-density  case  with  A/ 
=  50  and  Q  =  6.  In  these  computations,  the  temperature  was 
allowed  to  vary  while  the  density  was  maintained  constant. 
The  left-hand  side  of  the  figure  shows  the  results  at  the  e^ly 
stages,  where  the  fundamental  mode  dominates  the  dynamics, 
and  the  right-hand  side  shows  the  results  in  the  downstream 
part  where  pairing  is  the  dominant  mechanism  in  the 
dynamics.  The  results  in  this  figure  will  be  discussed  in 
connection  with  the  results  of  the  variable-density  cases. 

In  Fig.  7,  we  focus  on  the  fundamental  mode  and 
show  the  above  variables  in  the  section  1.43  <  x  <  2.43  at 
time  t  =  20  foT  Af  =  50  and  Q  =  4  and  6.  In  Fig.  8,  we 
concentrate  on  the  subharmonic  mode  and  show  the  same 
variables  in  the  section  2.82  <x<  3.82  at  time  r  =  15.5  for  Af 
=  50, 100  and  Q  =  4.  In  general  we  find  that  increasing  the 
amount  of  heat  release  by  using  larger  values  of  Q  has  a 
stronger  effect  on  the  dynamics  of  the  layer  than  increasing  the 
rate  of  energy  release  which  is  achieved  by  using  larger  Af. 
This  is  exhibited  by  the  extra  flattening  of  the  eddies,  and  then- 
further  expansion  in  the  streamwise  direction,  as  the  heat 
release  increases  by  increasing  Q  in  Fig.  7  than  by  using 
higher  values  of  Af  in  Fig.  8.  In  both  figures,  increasing  the 
rate  of  heat  release  leads  to:  (a)  flattening  the  eddies  although 
the  total  cross  section  remains  the  same  or  increases;  (b) 
decreasing  the  size  of  the  out-of-phase  eddy  with  respect  to 
that  of  the  in-phase  eddy;  (c)  maintaining  higher  rates  of  heat 
release  in  the  braids  due  to  the  weaker  roll-up;  and,  (4) 
modifying  the  vorticity  distribution  inside  the  eddies 
extensively.  Both  the  pr^uct  mass-fraction  and  the  reaction 
rate  remain  high  within  the  braids  as  the  rate  of  energy 
deposition  increases  indicating  that  the  strain  there  is  reduced, 
thus  the  flame  is  not  extinguished  to  the  same  degree  as  in  the 
uniform-density  case.  In  the  second  section,  the  fate  of  the 
out-of-phase  eddies,  as  they  are  strained  and  gulped  by  the 
larger,  in-phase  eddies  is  also  shown.  The  relative  location  of 
the  eddies  is  much  different  than  it  was  in  the  uniform  density 
case. 

The  similarity  between  the  product  mass-fraction  and 
the  vorticity  is  observed  in  these  results  albeit  to  a  lesser 
degree  than  what  we  observed  in  the  uniform-density  case. 
This  is  because  vorticity  generation  due  to  the  baroclinic 
torque,  produced  as  a  thin  zone  of  low-density  fluid  forms 
within  the  reaction  zone,  substantially  distorts  the  original 
vorticity  within  the  layer.  The  area  where  the  vorticity  is  most 
disturb^  by  the  combustion  process  is  where  the  reaction 
zone  is.  There,  the  vorticity  is  intensified  around  the  reaction 
zone  on  one  side  of  the  ^dy  and  changes  sign  across  the 
teaction  zone  on  the  other  side  of  the  eddy.  This  indicates  that 
baroclinic  vorticity  plays  a  significant  role  in  the  developing 
structures  of  the  shear  layer.  The  complex  structure  of  the 
vorticity  field  resembles  qualitatively  that  observed  in  the  case 
of  a  nonreacting,  variable-density  shear  layer.  It  is  difficult  at 
this  stage  to  distinguish  between  the  role  played  by  the 
baroclinic  vorticity  and  volumetric  expansion  since  they  are 
strongly  connected  with  the  reaction. 

The  plots  of  the  reaction  rate  show  that  the  teaction 
zone  becomes  almost  a  ring  around  the  large  eddies.  While  it 
is  clear  that  at  the  early  stages  the  reaction  starts  inside  the 
eddies  where  the  mixed  fluid  exists,  the  reaction  zone  moves 
outwards  in  the  large  eddies  as  the  original  mixed  reactants  are 
burned  and  are  replaced  by  two  newly  entrained  streams.  The 
same  reaction-zone  structure  is  observed  in  the  downstream 
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part  of  the  domain  where  the  interaction  between  the 
fundamental  eddies  produce  larger  structures.  Clearly,  heat 
release  does  not  affect  the  relative  position  of  the  reaction  zone 
to  the  vorticity  structures  as  it  changes  the  latter. 


V.  CONCLUSIONS 

Reacting  shear  layer  simulations  using  the  transport 
element  method  has  been  used  to  investigate  the  effect  of  heat 
release  on  the  flow  dynamics  in  a  spatially  developing  flow  of 
nonpremixed  reactants.  Results  were  obtained  for  a  uniform- 
density  reacting  shear  layer  and  three  cases  of  variable-density 
reacting  shear  layers  with  different  enthalpy  of  reaction  and 
reaction  frequency  factor.  In  all  cases,  a  single  step, 
Arrhenius  chemical  reaction  was  used  to  model  the  combustion 
process.  In  all  of  our  calculations,  forcing  was  used  to  initiate 
the  roll-up  process.  This  is  because  we  found  that  without 
forcing,  heat  release  may  completely  suppress  the  instability 
leading  to  roll-up. 

Results  indicate  that  the  overall  growth  rate  of  the  shear 
layer  is  reduced  due  to  the  heat  release  when  density  variation 
is  considered.  Increasing  the  amount  of  energy  deposited  is 
more  effective  in  reducing  the  growth  rate  than  the  rate  of 
energy  deposition.  In  the  reacting  shear  layer,  as  the  rate  of 
energy  deposition  increases,  the  rate  of  growth  and  the  rate  of 
product  formation  increase.  The  size  of  the  eddies  increases 
due  to  the  volumetric  expansion,  and  the  strain  field  induced  in 
their  braids  becomes  smaller  than  that  experienced  in  the 
uniform  density  case.  Thus,  while  the  reaction  rate  inside  the 
eddies  is  reduced  due  to  the  slower  roll-up,  it  is  increased  due 
to  the  reduction  of  strain.  Pairing,  a  process  responsible  for 
the  growth  of  the  shear  layer  downstream  the  roll-up  stage,  is 
suppressed.  This  is  due  to  the  change  of  phase  relationship 
among  neighboring  eddies  due  to  the  acceleration  associated 
with  heat  release.  Pairing  is  replaced  by  a  process  of 
tearing/gulping  in  which  in-phase  eddies  exert  a  strong  strain 
on  out-of-phase  eddies  and  absorb  their  vorticity. 
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Figure  2.  The  forcing  function  used  to  accelerate  the  initiation 
of  roll-up  and  pairing.  The  displacement  shown  in  the  figure 
represents  the  location  of  the  computational  elements 
immediately  downstream  the  splitter  plate. 


7 


Ar-  !0.0  Qa-  6.0  (•  29  00  ODCKTS*  4773  Ar-  S0  0  Oo-  6  »  f  20  «0  EIBCNTS*  MCS 


Figure  3.  The  location  and  velocity  vectors  of  the  voitex  elements  of  the  uniform-density  (left)  and 
variable-density  (right)  reacting  shear  layers,  both  with  Af=  SO  and  Q  =  6. 
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Figure  5.  The  growth  of  the  reacting  shear  layer  shown  in  terms  of  the  product  thickness  (left)  and 
the  reaction  thickness  (right). 


Figure  6.  From  top,  the  product  mass-fraction,  reaction  rate  and  vonicity  for  the  reacting  shear 
layer  at  /!/=  50  and  Q  =  6  for  the  uniform  density  case,  shown  for  1.43  <  x  <  2.43  (left)  and  2.45 
<  X  <  3.45  (right).  In  all  cases  0.25  <  y  <  0.75. 


Figure  8.  from  top,  the  product  mass-fniction  and  reaction  rate  fotQ  =  4  and  Af  =  50  (left)  and  Af 
=  100  (right).  The  results  are  shown  for  2.82  <  x  <  3.82.  In  all  cases  0.25  <  y  <  0.75. 
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